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FOREWORD 

This  report  was  prepared  by  Mr.  Windcll  F.  Ingrain  of  the  Systems 
Software  Section,  Operations  Branch,  Automatic  Data  Processing  Center, 

U.  S.  Army  Engineer  Waterways  Experiment  Station.  The  report  is  es¬ 
sentially  a  thesis  submitted  by  Mr.  Ingram  in  partial  fulfillment  of 
the  requirements  for  the  degree  of  Master  of  Science  in  Engineering  to 
the  Faculty  of  Mississippi  State  University,  and  is  a  study  concerned 
with  vehicle  dynamics.  The  study  described  herein  was  conducted  under 
the  auspices  of  the  Mobility  and  Environmental  Systems  Laboratory,  Water 
ways  Experiment  Station  under  DA  Project  lTl62112A0i*6,  "Traf ficability 
and  Mobility  Research,"  Task  03,  "Mobility  Fundamentals  and  Model 
Studies,"  under  the  sponsorship  and  guidance  of  the  Research,  Develop¬ 
ment  and  Engineering  Directorate,  U.  S.  Army  Materiel  Command.  The 
st,udy  was  accomplished  under  the  general  direction  of  Messrs.  S.  J. 
Knight  and  V/.  G.  Shockley. 
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publication  of  this  report.  Mr.  F.  R.  Brown  was  Technical  Director. 
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CHAPTER  I :  INTRODUCTION 


A  vehicle  body  considered  as  rigid  possesses  six  degrees  of 
freedom,  and  may  experience  any  of  six  different  types  of  motion. 
Motion  in  the  direction  of  travel  of  the  vehicle  is  related  to  accel¬ 
eration  and  braking  characteristics  and  is  generally  called  the 
'performance"  motion.  Yawing  and  sideslipping  are  associated  with 
steering  control  and  are  called  "handling"  motions.  Bounce,  pitch, 
and  roll  are  considered  to  be  "ride"  motions  of  the  vehicle. 

In  a  general  sense,  the  analysis  of  ground  vehicle  ride  can  be 
viewed  as  a  problem  involving  the  riding  surface,  the  vehicle,  and 
the  vehicle  occupants.  A  complete  analysis  of  occupant  comfort 
would  deal  not  only  with  the  physical  phenomena  stimulating  occu¬ 
pants,  but  also  with  the  physiological  and  psychological  responses 
of  the  occupants.  The  vehicle  occupants  are  stimulated  by  vehicle 
notions  and  by  noises  and  vibrations  coming  from  the  power  plant, 
the  power  transmission  system,  and  other  sources.  While  stimulated 
by  the  vehicle's  response  to  steering,  acceleration,  and  braking,  the 
occupants'  most  frequent  stimulation  is  due  to  the  notions  of 
bounce,  pitch,  and  roll.  These  components  of  ride  dynamics,  which 
may  be  defined  as  the  vehicle's  response  motions  to  excitations 
from  riding  surface  irregularities,  are  the  only  vehicle  motions 
considered  in  this  study.  While  some  researchers  consider  roll 
motions  to  lie  in  the  domain  of  handling  rather  than  ride  charac¬ 
teristics,  roll  has  been  shown  to  contribute  significantly  to 
vertical  motions  of  the  vehicle  occupants,  especially  in  cross-country 
vehicles  where  terrain  profiles  along  the  paths  of  the  two  sides  of 


the  vehicle  can  be  significantly  different.  The  effects  of  engine 


vibrations  and  vehicle  frame  flexure  are  not  generally  considered  to 
be  significant  contributors  to  ride  dynamics  and  will  not  be  con¬ 
sidered  here . 

Generally,  the  primary  objective  of  an  analysis  of  ride 
dynamics  is  to  improve  the  ride  characteristics  of  vehicles  through 
an  understanding  of  the  mechanics  of  vehicle  response  to  riding 
surface  irregularities.  However,  some  analyses  may  be  for  the 
specific  purpose  of  determining  which  existing  vehicle  design 
exhibits  most  favorable  response  characteristics  to  u  given  type 
riding  surface,  thus  facilitating  selection  of  a  vehicle  for  a 
particular  application. 

The  three  methods  used  in  ride  analysis  are  generally  referred 
to  as  observation,  experimental  analysis,  and  mathematical  analysis. 

Of  the  three  methods,  the  oldest  and  most  obvious  is  that  cf 
observation.  This  is  simply  to  observe  the  ride  of  the  vehicle  on 
various  types  of  riding  surfaces,  make  changes,  and  observe  their 
influence.  This  has  been  the  most  widely  used  method  and  it  was 
largely  through  the  use  of  this  method  that  the  ride  of  our 
present  vehicles  evolved.  The  observation  method  has  several 
drawbacks : 

1.  The  method  is  subjective,  depending  on  the  observer  for  memory 

and  evaluation,  and  does  not  give  a  numerical  measurement  of 
ride  characteristics. 

2.  While  detecting  disturbing  ride  characteristics,  the  observer 

has  no  way  to  determine  the  source  of  the  disturbance  or 
what  changes  are  necessary  to  improve  the  condition. 

3.  Any  changes  in  vehicle  design  must  be  fully  implemented  in 

the  prototype  vehicle  before  evaluation  can  be  made. 


Experimental  methods  improve  on  the  observation  method  by 
employing  instrumentation  to  quantify  the  vehicle  response.  Modern 
electronic  instrumentation  makes  it  possible  to  measure  and  record 
many  vehicle  motions  plotted  against  time.  Experimental  methods, 
which  are  widely  used,  share  with  the  observation  method  the  dis¬ 
advantage  of  requiring  a  vehicle  prototype  for  design  evaluation. 
Proposed  design  changes  can  be  evaluated  only  by  construct  Lon  of 
the  changed  component  and  implementation  in  the  prototype.  It 
is  also  necessary  to  construct  a  test  track  or  in  some  manner 
simulate  the  riding  surface  for  which  the  vehicle  is  intended. 

Mathematical  analyses  have  been  made  feasible  by  the  advent  of 
high-speed  analog  end  digital  computers.  Potential  advantages  of 
mathematical  analyses  are 

1.  They  provide  a  fundamental  understanding  of  t:  rrc-lation 

of  the  various  components  of  the  vehicle.  This  ; reviles  a 
better  understanding  of  the  influence  of  a  design  change. 

2.  The  effect  of  a  proposed  design  change  can  be  es'.ima'ed 

without  building  a  prototype. 

3.  Comparison  between  many  existing  vehicle  designs  on  varying 

riding  surfaces  can  be  made  by  changing  the  parameters 
associated  with  the  mathematical  model. 

Experimental  analyses  are  often  used  in  conjunction  with  and  for 

verification  of  mathematical  models. 

Historical  Background 

Attempts  at  analytical  approaches  to  vehicle  dynamics  date  back 
to  shortly  after  construction  of  the  first  automobile.  One  of  the 
first  by  Rowell  (fl)*  in  1922  dealt  with  undamped,  free  vibrations. 


* Numbers  in  parentheses  refer  to  similarly  numbered  entries  in  the 
Bibliography. 
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In  1926  Guest  (9)  developed  a  system  of  differential  equations  to 
describe  the  vibration  motion  of  a  vehicle.  His  equations  included 
effects  of  body  bounce,  body  pitch,  and  axle  bounce  in  a  system 


P 

k 


consisting  of  body  mass,  suspension  springs,  axle  masses,  and  tire 
vibration  for  body  bounce  and  pitch  to  the  motions  felt  by  passen¬ 
gers  located  in  different  parts  of  the  vehicle.  In  193'1  01  ley  (l8) 
extended  the  work  of  Guest  to  include  independent  suspension  systems. 
Through  calculations  and  experiments  he  determined  the  effects  of 
changes  to  a  vehicle's  design  parameters  (suspension  and  damping 
rates)  on  the  vehicle's  periods  of  free  vibration  in  body  bounce  and 
pitch.  He  also  determined  that  ride  comfort  is  not  solely  a  matter  of 
opinion,  but  thut  a  definite  determination  of  the  motion  functions 
that  produce  discomfort  could  be  made  with  variation  being  small 
between  individuals.  His  work  made  it  possible  to  visualize  clearly 
many  of  the  problems  associated  with  vehicle  springing.  In  193r 
James,  et  al.  (ll),  using  the  analysis  by  Guest  as  a  basis,  developed 
the  first  method  for  obtaining  the  pitch  moment  of  inertia  of  a 
vehicle  body  experimentally.  This  was  a  significant  contribution 
since  an  uccurate  value  of  the  moment  of  inertia  is  essential  to 
mathematical  analysis  of  ride.  These  early  attempts  at  ride  analysis, 
while  producing  considerable  insight  into  the  problem,  soon  became 
too  mathematically  complex  for  solution  by  manual  methods  and  it 
became  apparent  that  some  sort  of  simulation  of  the  vehicle  would  be 
necessary . 

In  19**3  Schilling  and  Fuchs  (23)  built  a  mechanical  differential 


analyzer  to  solve  the  nonlinear  differential  equation  associated  with 


a  single  mass,  single  spring  damper  system.  Their  analyzer  permitted 
the  inclusion  of  u  nonlinear  characteristic  for  the  shock  absorber 
and  was  used  to  determine  the  transient  vertical  accelerations 
that  might  be  experienced  by  a  vehicle  passenger.  They  demons i rated 
how  these  accelerations  could  be  altered  drastically  by  changes  in 
the  characteristics  of  the  shock  absorber.  This  analyzer  was  con¬ 
sidered  to  be  the  first  analog  computer  used  in  suspension  analysis 
and  was  a  significant  advance  to  the  state  of  the  art. 

In  the  early  1950's  there  was  considerable  work  directed 
toward  passenger  car  suspension  analysis  in  the  uutornotive  industry 
as  a  result  of  the  advent  of  *  he  analog  computer.  A  significant 
advance  was  made  in  1953  by  Jeska  (12)  with  the  presentation  of  a 
model  that  included  vertical  motions  of  tee  front  and  rear  wheels 
together  with  the  bounce  and  pitch  of  the  b >dy .  The  unique  contri¬ 
bution  here  was  that  the  analysis  included  at.  actual  road  profile, 
measured  by  a  photographic  technique,  as  a  dri  ring  function  to  the 
..  del . 

In  1955  Bodeau,  et  al.  (2),  presented  a  very  good  general 
i'seussion  of  ride  analysis  and  a  detailed  analysis  using  a  nine- 
.  gree-of- freedom  model  to  describe  a  passenger  car.  This  work, 
wh.ich  included  both  an  experimental  analysis  and  a  mathematical 
analysis  using  an  analog  computer,  attempted  to  determine  human 
response  characteristics  in  terms  of  passenger  comfort  limits.  Up 
to  this  point,  little  had  been  learned  regarding  passenger  tolerance 
to  vibrations. 

Kohr  (10)  in  I960  constructed  a  servo-controlled  motion  simu¬ 
lator  which  could  seat  two  occupants  and  which  followed  motions 


determined  by  an  analog  computer.  This  "ride  simulator"  included 
seven  degreesof  freedom,  i.e.,  body  bounce,  roll,  and  pitch,  plus 
bounce  and  roll  of  each  of  the  two  axles.  A  road  profile  was 
recorded  on  magnetic  tape  and  then  "played"  into  the  analog  computer 
cvS  input  disturbance  to  the  mode],  suspension  system.  The  simulator 
was  occupied  by  a  driver  and  passenger  to  provide  a  means  of  deter¬ 
mining  the  effects  of  the  vibrationaj.  characteristics  of  humans. 

The  ride  simulator  was  an  important  advancement  for  the  automotive 
.idustry  because  it  enabled  improvements  in  passenger  car  ride  to  be 
made  with  a  minimum  of  time  and  effort. 

Up  until  the  late  1950's  most  of  the  research  in  ride  dynamics 
had  been  directed  at  on-the-road  vehicles  such  as  passenger  cars. 
However,  in  this  case  most  of  the  severe  effects  of  ride  motions  had 
been  overcome  by  altering  the  environment  rather  than  the  vehicle; 
i.e.,  we  have  designed  and  built  smooth  roads  that  produce  rather 
minor  disturbances  to  the  suspension  system  as  compared  with  cross¬ 
country  locomotion. 

In  the  late  1950 's  the  military  began  to  recognise  the  signi¬ 
ficance  of  the  ride  dynamics  problem  in  off-road  mobility.  For 
off-road  locomotion  to  achieve  maximum  effectiveness,  ride  motions 
due  to  terrain  roughness  must  cause  neither  mechanical  failure 
nor  human  response  that  reduces  the  average  velocity  of  the 
vehicle.  Mechanical  failure  of  military  vehicles  due  to  terrain 
roughness  is  not  common,  thus  the  military  driver  in  an  emergency 
will  always  drive  to  the  limits  of  his  tolerance  to  vibratory 
motions.  Analysis  of  ride  stimuli  to  the  driver  and  crew  in  this 
environment  bears  little  resemblance  to  determination  of  comfort 


6 


levels  in  the  usual  automotive  sense.  The  objective  is  to  determine 
limiting  speeds  for  cross-country  vehicles  in  terms  of  the  vehicle's 
design  parameters,  the  terrain  roughness,  and  human  tolerance  to 
vibratory  motions. 

As  a  result  of  this  military  awareness  of  ride  dynamics, 
mathematical  models  to  assess  the  influence  of  various  design  al¬ 
ternatives  on  the  performance  of  military  vehicles  have  become 
important.  Greater  requirements  for  ground  mobility  have  produced 
requirements  for  new  vehicles.  There  is  no  longer  time  for  expensive 
"cut  and  try"  methods  in  search  of  new  designs.  Technology  has 
increased  the  number  cf  possible  engineering  choices  to  the  extent 
that  exploratory  building  and  testing  of  all  the  alternatives  is 
impossible.  This  has  led  to  the  recognition  that  mathematical  models 
of  terrain-vehicle  systems  are  necessary  to  explore  the  possible 
choices  and  select  a  rational  optimum. 

Thus  the  19b0's  saw  the  birth  of  several  studies  of  off-road 
vehicle  mobility,  including  not  only  analysis  of  ride  dynamics,  but 
also  analysis  of  trafficabili ty  and  maneuverability.  Van  Deusen  (PM 
conducted  a  comprehensive  study  of  the  ride  dynamics  aspect  of 
off-read  mobility  including  a  theoretical  treatment  of  the  dynamic 
behaviour  of  vehicles  and  an  analysis  of  human  response  to  vehicle 
vibration.  Pradko,  et  al.  (20),  conducted  basic  research  n  the 
effects  of  venicle  vibrations  on  man.  On  the  assumption  that  ride 
comfort  is  strongly  influenced  by  vertical  acceleration,  he 
developed  quantities  to  describe  human  dynamic  tolerances.  In  1965 
the  FMC  Corporation  (7)  conducted  a  study  for  the  Waterways  Experiment 
Station  (WES)  to  determine  the  feasibility  of  using  a  digital  computer 


i 

to  evaluate  the  ride  dynamics  of  a  vehicle  traversing  a  hard-surface 
terrain.  It  was  this  FMC  study  that  provided  a  starting  point  for 
the  author's  work  on  this  thesis.  In  1966  Cornell  Aeronautical 
Laboratory  (5)  began  a  comprehensive  study  of  ground  mobility,  the 
result  of  which  included  a  sophisticated  vehicle  dynamics  model  using 
a  digital  computer.  The  model,  however,  was  not  verified  against 
experimental  results.  This  reference  contains  a  comprehensive  list 
of  all  ground  mobility  studies  conducted  by  and  for  the  U.  S.  Govern¬ 
ment  in  recent  years.  The  efforts  in  off-road  mobility  research, 
both  completed  and  in  progress,  are  extensive. 

Two  basic  approaches  have  been  taken  in  attempting  to  interpret 
ride  vibrations  and  correlate  them  with  passenger  response  and 
terrain  characteristics.  The  deterministic  approach  uses  actual 
terrain  or  obstacle  displacement  records  as  a  function  of  time  as 
input  to  the  model  and  produces  information  concerning  vehicle  motion 
at  each  instant  of  time.  This  allows  determination  of  maximum  ac¬ 
celerations  and  displacements,  etc.  The  stochastic  approach 
involves  describing  the  terrain  profile  statistically  and  performing 
a  statistical  analysis  of  the  response.  The  output  of  a  stochastic 
analysis  consists  of  statistical  quantities  associated  with  the 
vehicle  motions.  Although  more  investigations  have  been  made  using 
the  deterministic  approach,  emphasis  has  apparently  shifted  to 
stocnastic  techniques. 

Virtually  all  models  developed  use  either  an  analog  computer 
or  a  digital  computer  in  the  analysis.  The  determination  of  which 
is  used  is  related  to  the  nature  of  the  study  and  the  computing 


equipment  available  to  the  researcher.  Generally,  the  analog  approach 


i 

r 


is  more  limited  in  number  of  degrees  of  freedom  than  the  digital 
approach  since  many  organizations  interested  in  vehicle  ride  dynamics 
do  not  have  a  large  scale  analog  computer  capable  of  including 
all  the  degrees  of  freedom  desired  in  a  complex  vehicle  simulation. 
Digital  computers,  though  readily  available  in  sizes  capable  of 
solving  many  degree  of  freedom  vehicle  simulations,  do  not  provide 
the  speed  of  solution  as  does  the  analog  and  often  involve  long 
program  run-times  for  representative  terrain  lengths. 

Purpose 

The  general  purpose  of  this  study  was  to  provide  an  improved 
capability  for  simulating  vehicle  ride  motions  felt  by  the  driver. 

The  specific  objectives  encompassed  by  the  study  were  essentially 
three.  First  was  the  development  of  a  practicable,  efficient 
digital  computer  simulation  model  of  the  three-dimensional  ride 
dynamics  of  a  single-hull  vehicle  with  an  arbitrary  number  of  axles. 

A  background  for  this  effort  was  provided  by  computerized  models 
developed  by  FI-'.C  Corporation  and  by  the  V/t’S.  The  second  objective 
was  validation  of  the  model  by  comparison  of  model  responses  to 
field  test  measurements  for  a  representative  vehicle.  The  final 
objective  was  to  investigate  the  effects  of  numerical  ir.tegratior 
step-size  on  the  model  responses  and  to  formulate  Titeria  for 
dynamic  managemen  of  the  step-size  to  achieve  optimum  model 
performance . 
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The  effort  consisted  entirely  of  developing,  implementing,  and 
evaluating  a  digital  computer  simulation.  Several  field  tests  for  a 
representative  vehicle  were  simulated  and  comparisons  are  presented 
in  t-he  form  of  graphs  of  meusured  and  simulated  vehicJe  responses. 
The  effects  of  integration  step-si zo  are  presented  in  the  form  of 
graphs  comparing  simulated  vehicle  responses  using  various  inte¬ 
gration  intervals.  Criteria  for  management  of  tne  integration 
interval  are  established. 

The  scope  of  the  study  did  not  include  the  execution  of  field 
tests  or  the  measurement  of  vehicle  parameters.  The  field  test  data 
were  provided  by  the  Mobility  Research  branch  (MRB)  of  the  Wr.C  and 
the  vehicle  parameters  used  in  the  simulation  came  from  the  MRB  and 
the  FMC  Corporation.  The  iatu  were  accented  as  being  as  reliable 
as  any  that  were  availabie  during  ho  course  of  this  study. 
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CHAPTER  II:  THE  MATHEMATICAL  MODEL 


The  Differential  Equations  of  Motion 

The  ride  motions  of  a  generalized  three-dimensional  vehicle  tra¬ 
versing  an  unyielding  riding  surface  can  be  described  by  writing  a  set 
of  simultaneous  differential  equations  relating  the  accelerations  of 
the  vehicle  components  to  the  spring  and  damping  forces  applied  by 
the  pneumatic  tires  and  suspensions  (springs  and  shock  absorbers). 
Considering  the  vehicle  body  and  axles  as  rigid  masses,  the  motions 
to  be  included  are  body  bounce,  body  pitch,  body  roll,  axle  bounce, 
and  axle  roll.  No  fore  and  aft  or  transverse  vibratory  motions  are 
to  be  considered  here.  By  viewing  the  vehicle  body  and  axles  as 
displaced  from  their  static  equilibrium  position  (Figure  l)*,  we  can 
identify  the  spring  and  viscous  damper  forces  contributing  to  body 
bounce  as  shown  below: 

Force  produced  by  vertical  displacement  of  the  body: 

2n 

-X  (1) 

J=1 

Force  produced  by  vertical  motion  of  the  body: 

.  2n 

-X  l's  (2) 

J=1 

Force  produced  by  body  pitch  displacement: 


’‘Each  tire  is  depicted  in  Figure  1  as  a  single  spring-damper  assembly 
for  convenience  only.  As  actually  modeled,  stiffness  of  each  tire  is 
represented  by  many  radially  spaced  springs  as  will  be  shown  on  the 
following  pages . 
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Figure  1.  Schematic  of  the  Vehicle  Model 
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Force  produced  by  body  pitch  motion: 
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Force  produced  by  body  roll  displacement: 
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Force  produced  by  body  roll  motion: 
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Force  produced  by  vertical  displacement  of  the  axles : 
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Force  produced  by  vertical  motion  of  the  axles: 
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Force  produced  by  axle  roll  displacement: 
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Force  produced  by  axle  roll  motion: 
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Summing  suspension  forces  and  equating  to  the  body  accelerating 


force  yields: 
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Similarly  the  differential  equation  describing  body  pitch 


can  be  written  by  summing  the  pitching  moments  about  the  body 


center 

of 

gravity. 
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Coupled  with  the  pitch  and  bounce  equations  is  an  equation 
describing  the  rolling  motion  of  the  vehicle  body.  This  is  obtained 
by  summing  the  rolling  moments  about  the  body  center  of  gravity. 
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Each  axle  assembly  requires  two  diffeiential  equations, 
one  describing  axle  bounce  and  one  describing  axle  roll,  to 
couple  with  the  three  body  equations.  The  equation  describing 


axle  bounce  is  obtained  by  summing  forces  from  the  suspensions 
and  tires  connected  to  an  axle.  The  axle  bounce  equation  is 
written  as: 


(14) 
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The  general  equation  for  axle  roll  is  obtained  by  summing 
the  rolling  moments  about  the  axle  center  of  gravity. 
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The  coupled  set  of  second  order  differential  equations 
describes  all  the  components  of  vehicle  motion  associated  with 
ride  dynamics,  i.c.,  body  bounce,  body  pitch,  body  roll,  axle 
bounce,  and  axle  roll,.  Thus  for  un  n-axled  vehicle  the  system 
described  has  N  degree  of  freedom  where 

N  =  2n  +  3 

It  includes  all  the  nonlinearities  associated  with  large 
rotational  motions,  i,e„,  sinO,  cos9,  s i n4> ,  cos$,  etc. 

The  Segmented  Tire  Concent 

Many  vehicle  dynamics  models  use  a  representation  of  the 
pneumatic  tire-riding  surface  contact  mechar. ism  similar  to  that 
shown  in  Figure  2.  The  flexure  property  of  the  tire  is  represented 
by  a  single  spring  element  which  may  have  either  linear  or  non¬ 
linear  properties.  Dissipation  of  energy  is  accommodated  by  u 
viscous  damper  which  may  also  be  linear  or  nonlinear.  Using  this 
technique,  the  single  point  of  contact  is  deflected  by  changes  in 


elevation  along  the  riding  surface  profile. 


Figure  2.  Depiction  of  a  Point-Contact  Concept 


This  point-contact  method  may  be  adequate  in  many  cases 
where  the  riding  surface  is  relatively  smooth  and  there  are  no 
abrupt  changes  in  surface  slope.  For  on-the-road  vehicle  models 
this  method  may  be  sufficiently  representative  of  the  contact 
mechanism  while  providing  a  relatively  simple  implementation. 

For  off-roud  inodeis  where  large  surface  irregu]  arities , 
obstacles,  and  abrupt  slope  changes  may  be  encountered,  the  point- 
contact  concept  is  a  questionable  representation  of  the  tire- 
terrain  interface.  It  provides  none  of  the  "enveloping"  charac¬ 
teristics,  .  f  a  pneumatic  tire.  That  is,  small  terrain  features 
cunr.ot  produce  deflections  of  small  areas  of  the  tire  carcass. 
Every  deflection  of  the  point  of  contact  produces  forces  repre¬ 
sentative  of  a  deflection  of  the  entire  contact  surface  by  the 
same  magnitude.  Figure  i  snows  an  example  of  this  condition. 

The  small  round  obstacle  produces  a  deflection  at  the  contact 
point  that  produces  forces  of  the  same  magnitude  as  would  be 
produced  by  a  deflection  of  the  entire  contact  surface  as  shown  by 
the  solid  line.  .'he  obstacle  would  actually  >e  somewhat  enveloped 
by  tile  tire,  the  tire  surface  being  more  ciosely  approximated  by 
t.ue  dashed  line. 

Seeking  a  more  realistic  tire  simulation  method,  A.  S.  Lessem 
•  1  i )  of  the  VIES  Mobility  Research  Branch  conceived  a  variable-area 
contact  mechanism  that  simulates  tire  envelopment  properties  and 
provides  for  transmission  of  horizontal  as  well  as  vertical  forces 
through  the  tire.  The  flexure  property  of  the  tire  is  represented 
by  many  radial  springs  spaced  throughout  an  "active  region"  of 
the  tine,  i.e.,  that  area  of  the  tire  that  may  contact  the  riding 
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Figure  3.  Example  of  Excess  Deflection  Produced  by  Point-Contact 
Method 
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surface  (Figure  1 ) .  Each  spring,  representing  one  "segment"  of 
the  tire  carcass,  was  considered  linear  with  spring  coefficient  K 
being  the  some  for  each  segment.  Individual  springs  of  the  segmented 
tire  are  deflected  by  the  riding  surface  profile;  each  spring 
deflection  produces  o  force  vector;  summing  these  force  vectors 
produces  a  resultant  force  vector  representative  of  the  spring 
force  transmitted  to  the  vehicle  axle.  As  forward  vehicle  motion 
is  simulated  by  moving  the  nonrotating  active  region  relative  to 
the  riding  surface  profile,  surface  irregularities  Droduce  spring 
deflections  in  sequence  from  the  leading  oour.dury  to  tiie  trailing 
boundary  of  the  active  region  (Figure  5). 

Tc  implement  mis  method  of  simulating  tire  deflection,  it  was 
necessary  to  first  establish  a  segment  spring  constant  (K)  for  each 
tire  at  each  pressure  under  consideration.  An  experimentally  de¬ 
termined  static  force-deflection  relationship  for  the  tire  was 
obtained.  Then  u  single  point  on  the  force  deflection  curve  was 
chosen  for  the  purpose  of  computing  K.  It  was  found  that  in  most 
cases  a  vert ical  deflection  of  approximately  one  inch  was  an  optimum 
for  this  purpose.  From  Figure  It  static  equilibrium  requires  that 


n 

F  =2  Tka.cos$. 

“  i  1 

i=l 


Where : 

F  =  Vertical  load  at  the  chosen  deflection,  lb; 

K  =  Segment  spring  coefficient,  lb/in; 

=  Effective  radial  deflection,  i'*'*1  segment,  in; 

<J>  =  Force  direction  angle,  i^'1  segment,  radians; 
n  =  Number  of  segments  each  side  of  tire  center  line. 
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Then  K  can  be  computed  as 
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2  )  A. cos d>. 
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Figure  6  shows  an  example  of  the  actual  measured  force-deflection 


curve  for  a  specific  tire  and  the  force-deflection  curve  resulting 


from  a  segmented  tire  representation  using  twelve  segments  spaced  at 


ten-degree  intervals.  Note  that  the  modeled  deflection  curve  is 


piece-wise  linear  with  a  point  of  discontinuity  where,  as  deflection 


increases,  additional  segments  are  included  in  the  resisting  force. 


This  concept  was  first  modeled  dynamically  using  an  analog 


computer  to  simulate  the  motion  of  a  single  military  vehicle  tire. 


The  validity  of  the  results  was  verified  by  experimental  test  at 


the  WES. 


In  1968  this  mechanism  for  simulating  the  tire-riding  surface 


interface  was  incorporated  into  the  FMC  vehicle  dynamics  digital 


computer  program  by  J.  F.  Smith  (l6)  of  the  VIES  Automatic  Data 


Processing  Center.  Smith's  digital  implementation  of  the  segmented 


tire  concept  is  essentially  the  tire  mechanism  employed  herein. 


The  writer  has  made  modifications  to  the  tire  subprogram  for  the 


purpose  of  eliminating  some  oversights  in  the  original  code  and  to 


increase  the  speed  of  computation. 


Incorporation  of  the  segmented  tire  concept  into  the  general 


vehicle  dynamics  model  necessitates  modification  of  the  differential 


equations  describing  axle  bounce  and  roll.  In  equation  l4,  the  term 
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Figure  6.  Force  Deflection  Relationship  for  9.00-16  Tire  §  1*5  psi 
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is  replaced  by 


21 

j=2i-l 


(19) 


and  in  equation  15,  the  term 
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is  replaced  by 


KkkjLj 

j=2i-l 

2i 

1  V'j 

J=2i-1 


(20) 

(21) 


where  F  is  the  vertical  component  of  the  spring  force  vector 
J 

til 

transmitted  by  the  j  wheel  to  its  axle.  F  is  computed  by 

J 

summing  the  force  vectors  contributed  by  each  segment's  spring 
deflection.  The  spring  deflections  are  determined  by  computing 
the  point  of  intersection  of  each  segment  spring  with  the  riding 
surface  contacting  the  tire. 

Note  again  that  the  force-deflection  characteristics  of  the 
segmented  tire  model  include  many  discontinuities  in  stiffness  rates 
as  the  tire  carcass  flexes.  Effects  of  this  will  be  discussed 
later . 

While  no  attempt  has  been  made  to  incorporate  damping  properties 
into  each  tire  segment  as  with  spring  stiffness,  tire  damping  is 
included  in  the  digital  model  as  a  single  viscous  damper  responding 
to  motion  of  the  wheel  hub  (as  shown  in  Figure  2  with  the  point- 
contact  method). 
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Representation  of  the  Riding  Surface  Profile 

The  hard-surface  profile  which  serves  as  a  displacement 
function  for  the  tire  spring,  thus  in  essence  a  forcing  function  to 
the  system  of  masses,  springs,  and  dampers  representing  the  vehicle, 
is  represented  by  a  series  of  coordinates  connected  by  straight  line 
segments.  Elevations  along  the  surface  are  determined  by  simple 
linear  interpolation  between  the  points.  Examples  of  a  typical 
obstacle  test  profile  and  a  typical  segment  of  cross-country  terrain 
are  shown  in  Figure  7- 

One  might  consider  linear  interpolation  between  profile  points 
to  be  a  crude  method  of  describing  a  terrain  profile  for  the  purpose 
of  simulating  vehicle  traversal.  A  part  of  the  FMC  study  was 
directed  at  determining  practical  methods  for  describing  such 
profiles.  A  Fourier  series  was  used  as  one  method  of  mapping  a 
terrain  profile,  and  the  method  was  found  to  be  accurate  but  very 
inefficient  due  to  the  increase  in  computer  time  required  to 
generate  and  process  the  Fourier  series.  It  was  found  that  use 
of  linear  interpolation  between  profile  points  produced  sufficient 
accuracy  to  make  this  a  valid  approach  provided  that  the  profile 
point  spacing  was  not  excessive.  The  linear  interpolation  method 
represents  the  absolute  minimum  in  computer  time  that  would  be 
required  by  any  method. 


CHAPTER  III:  SOLUTION  OF  THE  EQUATIONS  OF  MOTION 


n 


Digital  simulation  of  a  dynamic  system  such  as  the  nonlinear 
system  described  by  the  equations  of  motion  for  the  generalized 
solid  axle  vehicle  requires  a  discrete  representation  of  the  con¬ 
tinuous  system,  i.e.,  a  discrete  state  model.  Given  an  initial 
system  state,  the  next  system  state  (after  some  time  interval  h) 
is  to  be  determined.  This  makes  the  system  simulation  a  problem 
in  numerical  integration. 

Numerical  integration  of  systems  of  second-order  differential 
equations  of  this  type  is  commonly  handled  by  first  expressing  the 
second-order  equations  as  a  system  of  first-order  equations  to 
which  general  numerical  integration  techniques  are  applicable. 
Letting  represent  a  general  displacement  of  a  vehicle  component, 
the  equations  of  motion  can  be  expressed  in  general  terms  as: 


yi  =  fi(t»yi>yi>y2’y2’--"’yr 


,yn)  ,i  =  l. . .n 


(22) 


"t  h 

where  i  denotes  the  i  differential  equation  and  n  is  the  number 
of  degrees  of  freedom,  i.e.,  the  number  of  differential  equations 
in  the  system.  Conversion  to  a  system  of  first-order  equations 
is  accomplished  by  substituting  z ^  for  y^ .  Thus  the  equations  of 
motion  can  be  expressed  as: 


y.  =  z±  (23) 

Zi  =  fi(t>y1>z1>y2’Z2,,0“yn’Zn)’1=1,’,jn 
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In  practice  this  means  that  computed  accelerations  (z  )  are 
used  to  solve  for  velocities  (z^,)  which  are  in  turn  integrated 
numerically  to  produce  displacements  (y^).  Given  the  initial 
values  of  displacement  and  velocity,  this  is  a  practicable  means 
of  integrating  dynamics  equations. 

Given  the  system  of  first-order  equations,  numerical 
integration  is  usually  accomplished  by  one  of  two  basic 
approaches.  One  is  a  two-step  predictor-corrector  procedure 
requiring  the  use  of  the  present  system  state  plus  several 
previous  states.  The  second  approach  is  a  one-step  procedure 
which  requires  computation  of  several  intermediate  values,  but 
requires  use  of  the  present  state  only,  e.g.  the  Runge-Kutta 
methods . 

A  cursory  review  of  the  literature  might  lead  one  to 
believe  there  is  an  abundance  of  good  numerical  integration 
techniques  and  that  the  choice  of  method  might  be  somewhat 
arbitrary.  However,  selection  of  the  integration  technique 
plays  a  large  part  in  the  overall  system  simulation.  Desired 
attributes  of  the  selected  method  include  speed,  accuracy, 
stability,  and  convenience  of  use.  Maximum  speed  car.  be  obtained 
through  selecting  a  fast  method  and  using  a  large  time  step; 
acceptable  accuracy  may  require  slower,  more  sophisticated 
methods  using  a  much  smaller  time  step. 

Fourth-order  Runge-Kutta  methods  are  among  the  most 
popular  and  widely  used  methods  of  today.  These  methods 


generally  have  the  following  advantages: 

1.  They  a’  c-  self-starting,  i.e.,  only  one  known  system  state  is 

requi  v i  to  initiate  the  procedure. 

2.  They  cU  not  usually  exhibit  numerical  instability  using 

properly  chosen  step  sizes. 

3.  Good  accuracy  is  usually  achieved  using  properly  chosen  step 

sizes . 

1*.  Relatively  small  amounts  of  core  memory  are  required  for 
processing . 

5.  They  are  convenient  to  use,  being  particularly  amenable  to 
general  purpose  subroutines . 

Disadvantages  generally  include: 

1.  No  estimate  of  accuracy. 

2.  No  way  of  knowing  if  the  step  size  is  adequate  without 

expensive  reruns  at  various  step  sizes. 

3.  Requirement  for  four  derivative  evaluations  per  integration 

step  compared  with  only  two  using  tie  fourth-order  predictor- 
corrector  methods. 


The  Runge-Kutta-Merson  Method 

The  Runge-Kutta-Merson  (RKM)  method  attempts  to  negate  seme  of 
the  disadvantages  associated  with  Runge-Kutta  methods  in  general  by 
using  a  fifth  derivative  evaluation  to  provide  an  estimate  of  dis¬ 
cretization  error  to  oe  used  for  automatic  step-size  adjustment. 

In  general  terms  the  Merson  formulae  are: 


■■  ‘  =  yf  +  (pii  +  +  I'si )/-  +  0(h/) 


where 


pn  =  nfi{tJ’yi’y2,*,-’yn)  /3 


2i  =  hfiltJ  +  (h/3)’^  +  pn>y2 +  pi 


(pi.) 

(25) 

(26) 
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In 
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P3i  =  hf.UJ  +  (h/3),  yj  +  (P11/2)  +  (P;?;)/2),yJ  (27) 

♦  <V2)  +  (V2)>— *  (V2>t<p?„/2)1/-' 


P^J  =  hr.{tj  +  (h/2),  yj  +  (BPjj/8)  +  (9Ij3l/8), 

+  (3P12/8)  +  (9P32/8 y|j  +  (3Pln/8)  +  (9Pln/8)}/3 
P5i  =  hfi{tJ+  h.yJ  +(3P11/u)-(9P31/2)+6Plii,yi  +(  3P. .  /2 ) 


-  (9P32/2)  +  6P1(2 


,yJ  +  ( 3P,  /2)-(9P,  12)  +  6f.  }/3 
n  In  3n  «n 


O(h^)  indicates  that  the  truncation  error  is  of  the 
order  h"  as  h  +  0,  This  term  is  not  a  part  of  the  calculations, 
but  merely  indicates  that  this  is  a  fourth  order  method.  An 
approximation  of  the  error  of  integration  (truncation  error) 
is  provided  by 


=  (P 


li 


-  (9P3i/2) 


+  It 


rfci 


IP  /2))/5 


which  is  generally  used  as  a  basis  for  automatic  step-sice 
aujustmt-'d- .  In  the  equations  above,  t  is  the  independent 
variable,  the  subscript  i  refers  to  the  differential 
equation,  the  superscript  J  refers  to  the  previous  system  state 
just  calculated,  while  j  +  1  is  the  system  state  to  be  calculated 


and  n  is  the  number  of  equations  in  the  system. 


As  with  all  Runge-Kutta  formulae,  Mer son's  equations  were 
derived  by  using  Taylor  series  expansions  and  algebraic  manipulations 
to  match  coefficients  of  terms.  The  process  always  results  in  more 
unknowns  than  equations,  requiring  an  arbitrary  selection  of  one  or 
more  constants.  This  arbitrary  selection  of  constants  has  resulted 
in  the  many  variations  of  the  basic  Runge-Kutta  method.  Merson's 
unique  contribution  was  in  his  determination  of  an  approximation  of 
the  truncation  error  inherent  in  his  formulae. 

The  nKM  process  can  be  viewed  in  geometric  terms  as  follows. 

For  the  curve  associated  with  the  i^1  differential  equation,  compute 
the  slope  { ^ / ( h/ 3 ) }  at  the  point  (y^,t^)  and,  using  it,  advance 
one-third  step  forward  and  compute  the  slope  {P^/(h/3)}  there. 

Using  the  average  of  the  first  and  second  slopes  computed,  start 
again  at  (y^,t‘^)  and  advance  one-third  step  forward  and  again 
sample  the  slope  { /(h/3)} -  Then  starting  again  at  (y^,t^)  and 
using  a  weighted  average  of  the  first  and  third  slopes  computed, 
advance  one-half  step  forward  and  compute  the  slope  (Pj^/(h/3)} 
at  the  half  interval  point.  Finally  using  a  weighted  average  of 
the  first,  third,  and  fourth  slopes  computed,  advance  forward  a 
full  time  step  from  (y^,t^)  and  compute  the  slope  { ^/(h/3)} 
at  the  point  (y^f\  t^+^).  Having  computed  all  the  intermediate 
derivatives,  establish  the  point  (y||+\  tl^  +  1)  by  advancing  one  full 
time  step  from  (y^,t^)  using  a  weighted  average  of  the  slopes 
computed  at  the  beginning,  the  midpoint,  and  the  end  of  the  time 
interval.  The  error  estimate  e^  (equation  30 )  can  be  viewed 
geometrically  as  an  evaluation  of  the  variation  of  the  slopes  at 
the  intermediate  point,  i.e.,  the  greater  the  variation  in  slopes 


t 


over  a  time  interval,  the  greater  the  error  introduced  in  that 
time  step.  Hence,  it  seems  reasonable  from  a  geometric  point  of 
view  to  use  some  evaluation  of  intermediate  slope  differences 
such  as  equation  30  as  an  indicator  of  the  need  for  automatic 
interval  adjustment. 

Management  of  the  Numerical  Integration  Step-Size 

The  primary  requirements  for  acceptable  performance  of  Merson's 
or  any  other  numerical  integration  technique  are  speed,  accuracy,  and 
stability.  All  three  factors  are  directly  coupled  to  the  choice  of 
the  timi  increment  employed.  Since  most  of  the  computing  time  involved 
in  simulation  of  vehicle  dynamics  (or  any  other  dynamic  system)  is  used 
in  evaluating  the  derivatives,  the  computing  cost  is  directly  propor¬ 
tional  to  the  number  of  derivative  evaluations  required,  hence  inversely 
proportional  to  the  size  of  the  time-step  chosen,  e.g.,  halving  the 
step-size  doubles  the  computing  cost.  In  complex  systems  such  as 
is  being  simulated  here,  the  corcpu' ing  time  required  per  time-step 
precludes  choosing  a  "safe,"  i.e.  very  small,  constant  time-step  to 
assure  accuracy  and  stability.  Such  a  choice  would  make  the  computing 
costs  so  high  as  to  render  simulations  of  an  effective  time  duration 
impracticable.  Choosing  a  large  constant  time-step  to  achieve 
reasonable  computing  costs  can  result  in  unpredictable  performance 
or  unreliable  results.  Thus  proper  management  of  the  time-step-size, 
automatically  by  the  computer  program,  can  mean  the  difference  between 
an  unreliable,  unpredictable,  or  impracticable  model  and  a  workable 


model  producing  acceptable  accuracy  at  an  acceptable  cost. 


What  constitutes  "management"  of  the  integration  time-step? 

Ideally,  the  objective  is  to  produce  accurate  results  at  a  minimum 
cost.  This  can  only  be  achieved  if  at  each  discrete  point  in  time, 
the  simulation  moves  forward  from  the  current  system  state  to  the 
next  system  state  choosing  the  largest  time  increment  that  maintains 
stability  and  acceptable  accuracy.  Time  step  management  thus  becomes 
the  automatic  determination  or  recognition  by  the  computer  program  of 
1.  The  initial  time-step 

?.  System  states  requiring  time-step  reduction 
3.  The  minimum  allowable  time-step 
h .  System  states  allowing  time-step  increase 
5.  The  maximum  allowable  time-step. 

The  literature  offers  very  little  guidance  in  these  mutters.  Many 
users  of  general  purpose  numerical  integration  packages  seem  to  choose 
the  integration  parameters  rather  arbitrarily  on  a  "try  some  numbers  and 
see  what  happens"  basis.  Realistically,  criteria  for  each  decision  must 
be  based  on  the  characteristics  of  the  system  being  simulated.  The  con¬ 
siderations  involved  in  establishing  time-step  management  criteria  for 
the  ride  dynamics  modei  will  be  discussed  here,  with  the  bulk  of  the 
results  of  the  investigation  into  this  area  being  presented  in  Chapter  V. 
Initial  time-step 

line  might  think  that  the  choice  of  the  initial  time-step  would 
be  rather  arbitrary  if  the  integration  procedure  had  automatic  interval 
adjustment  built  in.  However,  interval  adjustment  procedures  do  not 
automatically  choose  "the"  correct  time-step  at  a  given  point;  they 
simply  determine  if  a  change  in  interval  size  is  justified  based  <■>« 
system  responses  at  each  discrete  system  state.  A  given  system  state 
may  not  yield  a  time-step  increase  even  though  a  larger  interval 


35 


1.  The  advance  of  the  system  from  the  previous  discrete  point 

produces  one  or  more  deflections  or  velocities  that  exceed 
the  limits  of  the  tables  describing  the  characteristics  of 
the  vehicle  components,  e.g.,  the  maximur- possib' e  deflection 
of  a  spring  has  been  exceeded. 

2.  The  advance  of  the  system  produces  an  impossible  or 

unrealistic  geometric  configuration,  e.g.,  the  center  cr  a 
wheel  is  below  tne  riding  surface  or  the  angular  displacement 
of  the  vehicle  body  is  beyond  stable  limits. 

3.  Integration  errors  reported  by  the  Merson  formula  are  beyond 

prescribed  limits. 

Either  of  the  first  two  conditions  might  be  considered  justifica¬ 
tion  for  termination  of  the  simulation  on  the  grounds  that  the  con¬ 
ditions  might  be  produced  by  invalid  modeling  techniques  or  incorrect 
vehicle  parameters.  (The  original  FMC  model  terminated  upon 
encountering  the  first  condition  and  provided  no  checks 
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for  file  second.)  However,  these  system  states  can  also  be  encountered 
simply  because  the  integration  step  taken  is  too  large,  i.e.,  the 
.rimulution  has  taken  a  step  forward  too  big  for  the  conditions  en¬ 
countered  in  that  step  (possibly  severe  terrain  irregularities),  thus  not 
allowing  the  model  to  respond  to  the  external  stimuli  in  increments  small 
enough  to  simulate  actual  vehicle  response. 

In  using  the  integration  error  reported  by  the  Merson  formula, 
normal  procedure  is  to  compare  the  reported  error  with  some  pre¬ 
selected  maximum.  If  the  reported  error  exceeds  the  maximum  allowable, 
the  interval  is  reduced  as  stated  previously.  Since  the  reported 
error  is  ideally  of  the  order  h  ,  halving  the  time-step  theoretically 
reduces  the  truncation  error  by  a  factor  of  2^  or  32. 

Applying  Merson’ s  error  formula  to  the  vehicle  model  is  not  a 
straightforward  task.  Simulation  of  a  two-axle  vehicle  (seven  degrees 
of  freedom)  requires  the  integration  of  fourteen  first-order  equations, 
thus  the  possible  requirement  of  monitoring  fourteen  different  com- 
puteu  error  values  for  each  time-step.  Motions  of  the  vehicle  com¬ 
ponents  are  of  widely  varying  magnitudes,  frequencies,  and  velocities. 
Also  the  presence  of  discontinuities  inherent  in  the  representation 
of  tlie  vehicle  suspension  parameters,  the  terrain  profile,  and  the 
tire  produce  severe  changes  in  the  vehicle  characteristics  (in 
essence,  changes  in  the  system  being  modeled)  during  one  time 
interval,  especially  with  larger  intervals.  Hence,  there  are  several 
questions  regarding  the  applicability  of  Merson 's  error  formula 
to  this  model. 
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1.  Does  error  computed  by  the  Merson  formula  accurately  reflect 

degradation  of  the  vehicle  responses  due  to  excessive  time- 
step  size? 

2.  Which  integration  error  should  be  monitored? 

3.  Should  the  error  be  evaluated  on  an  absolute  or  a  relative 

basis,  i,e.,  should  the  magnitude  of  the  variable  being 
computed  be  considered? 

h.  What  constitutes  an  error  threshold,  i.e.,  that  error  value 

that  requires  a  time-step  reduction  to  maintain  satisfactory 
model  performance? 

These  questions  will  be  discussed  in  Chapter  V. 

Minimum  allowable  time-step 

A  minimum  allowable  time-step  is  normally  specified  to  prevent 
the  time-step  from  being  reduced  to  the  point  where  roundoff  error 
begins  to  dominate  the  solution,  resulting  in  numerical  instability. 
Roundoff  error  in  numerical  integration  can  be  minimized  if  the  time- 
step  is  chosen  such  that  the  derivative  multiplier  (h/3  for  Merson' s 
formulae)  can  be  represented  exactly  by  the  computer  using  a  minimum 
number  *f  binary  digits,  preferably  one.  Using  Merson's  method  this 
means  choosing  h  such  that  li  =  3(^n)  where  n  is  an  integer.  The 
resulting  multiplier  is  2n  which  in  computer  real  number  storage 
requires  only  one  significant  binary  digit.  Always  changing  the  time- 
step  by  a  factor  of  two  (halving  r  doubling)  maintains  this  condition. 

Considering  the  computer  used  in  this  study  (Honeywell  Gh h 0 ) 
the  number  of  signifies'  digits  that  can  be  represented  is  such  that 
for  the  relative  magnitude  of  derivatives  and  time  intervals  involved, 
roundcf f  is  of  little  consequence.  Thus,  the  minimum  time-step  can  be 
based  on  considerations  other  than  roundoff.  Realistically ,  it 
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ca:.  be  based  on  minimum  forward  motion  of  the  vehicle  in  one  time- 


step.  It  can  be  assumed  that  once  the  time-step  has  been  reduced  to 
the  point  that  forward  motion  of  the  vehicle  is  less  than  some  very 
small  distance  (quantified  in  Chapter  V),  further  redaction  of  the 
time-step  will  not  improve  the  result.  If  further  reduction  is  in¬ 
dicated  for  any  of  the  reasons  previously  specified,  the  simulation 
should  be  terminated  and  the  cause  of  the  malfunction  investigated. 
Time-step  increase 

Assuming  that  the  simulation  began  using  the  maximum  allowable 

time-step  and  at  some  point  the  interval  was  reduced  for  one  of  the 

three  reasons  specified,  then  at  some  later  point  (possibly  after  a 

severe  disturbance  has  been  passed),  the  time  step  could  be  safely 

increased.  The  only  indication  available  as  to  the  feasibility  of 

increasing  the  interval  size  is  the  integration  error  computed  by 

Kersou's  formula.  Again,  essentially  the  same  questions  arise  as 

previrusly  encountered  with  the  time-step  increase 

i.  Do  very  small  errors  justify  a  time-step  increase? 

Which  integration  error  should  be  used  as  a  basis  for  an 
increase? 

?.  L'hould  the  error  be  evaluated  on  an  absolute  or  relative 

basis? 

h.  What  constitutes  the  lower  error  threshold? 

It  could  be  argued  that,  ideally,  for  a  truncation  error  of  the 

order  h'  establishment  of  the  maximum  error  allowed,  E  ,  would  also 
’  max 

establish  the  lower  error  threshold,  E  .  ,  at  E  .  =  E  / 32 .  However, 

min  min  max 

this  would  produce  an  oscillating  situation  in  which  the  step  size  was 

constantly  being  alternately  increased  and  reduced.  A  more  effective 

value  might  be  E  .  =  E  /50.  Realistically,  the  computed  error 

min  max  r 


terra  is  only  a  rough  approximation  of  the  error  involved  and  it 


remains  to  be  determined  how  it  can  be  applied  to  the  vehicle 
simulation . 

Maximum  allowable  time-step 

It  would  be  unrealistic  to  allow  the  time-step  to  grow  without 
establishing  an  upper  bound.  If  allowed  to  grow  without,  limit,  if 
is  conceivable  that  during  traversal  of  a  relatively  smooth  riding 
surface  on  which  very  little  vibratory  motion  is  occurring,  the  time- 
step  could  grow  to  the  point  at  which  forward  motion  in  one  step  is 
such  that  significant  terrain  irregularities  or  obstacles  are 
completely  "jumped  over"  without  being  "felt"  at  all  by  the  model. 

It  is  also  well  to  consider  natural  frequencies  of  vehicle  components 
when  establishing  a  maximum  time-step.  One  might  argue  that  it.  is 
necessary  to  execute  a  minimum  of  eight  or  ten  integrations  per 
period  of  the  highest  frequency  motion  being  modeled.  Thus  the 

maximum  time-step  might  be  established  as  a  function  of  the 
forward  velocity  of  the  vehicle  or  of  the  period  of  vibration  of 
some  vehicle  component.  Criteria  will  oe  established  for  this  it. 
Chapter  V . 


Description  of  the  Computer  Program 

The  computer  program  (Appendix)  is  written  in  FORTRAN  IV  for  a 
Honeywell  G4^0.  Conversion  to  any  other  digital  computer  with  a  FORTRAN 
capability  should  be  relatively  simple  to  accomplish.  The  only 
problem  that  might  be  encountered  in  such  a  conversion  is  the  introduc¬ 
tion  of  significant  roundoff  error  in  some  machines  with  less  than  ten 
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significant  decimal  digits  in  its  reTc-esentation  of  floating  point 
numbers . 

The  program  begins  by  reading  user-provided  parameters  describing 
the  vehicle.  These  include  the  weight  of  the  body,  the  body  moment  of 
inertia  in  pitch  and  roll,  the  position  of  each  axle  relative  to  the 
body  center  of  gravity,  the  weight  and  roll  moment  of  inertia  of  each 
ax l.o  assembly,  the  position  of  each  suspension  and  each  tire,  and 
the  spring  and  damper  characteristics  of  each  suspension  and  tire.  The 
springs  and  shock  absorbers  are  described  by  force-deflection  and 
force-rate  tables.  The  tires  are  described  by  specifying  the  number 
of  segments  to  be  used,  the  angle  to  each  segment  center  line,  and 
the  spring  constant  for  each  segment.  The  tire  damping  characteristics 
are  read  as  force-rate  tables,  respectively.  In  preparing  force- 
deflection  tables  for  the  suspension  springs,  the  effect  of  limited 
relative  displacement  between  the  axle  and  body,  i.e.,  the  axle 
assembly  hitting  the  bump  stop,  should  be  provided  for  by  a  very  stiff 
force-deflection  relationship  beyond  that  point  of  maximum  spring 
deflection . 

The  program  next  reads  the  riding  surface  profile  which  is 
described  simply  as  (x,y)  coordinates  relative  to  the  initial  position 
of  the  vehicle.  A  different  profile  can  be  entered  for  each  side  of 
the  vehicle  so  that  terrain  can  be  depicted  as  actually  encountered, 
which  may  be  considerably  different  for  the  right  and  left  tires. 

If  the  static  equilibrium  configuration  for  the  vehicle  is  known 
(the  deflection  of  ei  ch  tire  and  suspension  spring  at,  static 
equilibrium),  this  can  be  read  in  and  the  program  can  proceed  di¬ 
rectly  to  the  simulation  of  forward  notion.  If  the  stutic  equilibrium 
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configuration  for  the  vehicle  is  not  known  or  if  it  is  desired 
to  simulate  a  vertical  drop  of  the  vehicle  (with  no  forward  notion), 
then  an  arbitrary  initial  configuration  can  be  read  in  and  the 
program  will  provide  simulation  of  vehicle  oscillation  until  all 
motion  is  damped  out  and  static  equilibrium  is  reached.  This  is  a 
rather  time-consuming  procedure  and  should  be  avoided  once  the 
static  equilibrium  configuration  has  been  recorded  for  a  given  sot 
of  vehicle  parameters. 

Having  been  given  the  static  equilibrium  configuration,  the 
program  reads  the  initial  forward  velocity  and  the  forward  accel¬ 
eration  desired.  (This  latter  must  be  constant  and  is  usually 
zero.)  Next,  reading  the  stop  distance  or  the  stop  time  completes 
the  definition  of  the  simulation  required. 

Then,  simulation  of  vehicle  responses  to  forward  motion  over 
the  specified  riding  surface  begins  and  continues  until  stop  distance 
or  stop  time  is  reached.  The  procedure  essentially  consists  of 
determining  the  forces  on  the  vehicle  masses  produced  by  the  current 
displacements  and  velocities  (system  state),  relating  forces  and 
masses  to  determine  accelerations,  and  integrating  accelerations 
successively  t.o  determine  velocities  and  displacements  at  the 
next  point  in  time.  This  new  system  state  produces  new  forces  and 
the  procedure  begins  anew.  The  time  step  is  adjusted  if  required  as 
described  previously.  All  system  responses  at  each  point  in  time  are 
recorded  on  permanent  storage  devices  so  that  retrieval  for  generating 
graphical  displays  or  for  statistical  analyses  is  possible. 


42 


CHAPTER  IV:  VALIDATION  FOR  A  REPRESENTATIVE  k  X  1j  VEHICLE 

To  establish  the  validity  of  the  model,  it  was  decided  to 
simulate  a  military  vehicle  for  which  measured  vehicle  responses 
were  available  from  field  tests  previously  conducted  by  the  MRB  of 
the  WES,  Vicksburg,  Mississippi.  The  vehicle  chosen  was  the 
M37 — a  3/lj-ton,  four-wheel  drive  truck  used  to  transport  personnel 
or  light  cargo.  In  the  field  tests  the  truck  was  instrumented  with 
accelerometers  to  record  vertical  acce 1  orations  at  three  points — on 
the  front  axle  near  the  right  wheel,  at  the  body  center  of  gravity, 
and  on  the  truck  body  under  the  driver's  seat.  It  was  also  equipped 
with  a  gyro  to  record  the  pitch  motion  of  the  vehicle. 

Vehicle  Parameters 

The  vehicle  parameters  used  in  the  simulation,  with  the  exception 
of  the  tire  data,  were  taken  from  information  provided  to  WES  by  the 
FMC  Corporation.  It  is  known  that  the  force-deflection  relation¬ 
ships  for  the  springs  and  the  force-rate  relationships  for  the  shock 
absorbers  were  measured  from  disassembled  vehicle  components  and 
that  the  parameters  describing  the  relative  geometric  position  of 
the  vehicle  axles,  suspensions,  tires,  etc.,  were  measured.  It  is 
thought  that  the  moment  of  inertia  values  are  computed  values  and 
are  of  questionable  accuracy.  The  tire  segment  spring  stiffness 
was  computed  as  described  in  Chapter  II  using  measurements  made  by 
the  MRB.  Tire  damping  was  computed  by  the  MRB  based  on  drop  tests 
conducted  with  a  single  wheel.  The  tire  parameters  are  considered 


reliable . 


Courses  Traversed 

Four  field  tests  were  chosen  for  simulation-three  involving 
single  obstacle  traversal,  only  and  one  cross-country  test.  The 
obstacle  tests  were  for  a  single  obstacle  sine  and  shape — a  six- 
inch-high,  thirty-degree  ramp-type  wooden  obstacle.  The  obstacle 
data  used  represented  obstucle  impact  at  three  different  speeds — 
5.65  ft. /sec.,  15.7  ft. /sec.,  and  2P..0  ft. /sec.  The  cross-country 
run  simulated  consisted  of  traversal  of  300  feet  of  hard,  rocky 
terrain  located  at  a  WES  test  site  in  Nevada.  The  test  was  run  at 
approximately  eight  miles  per  hour. 

Predicted  and  Measured  Responses 

Figures  8,  9,  and  10  show  the  measured  and  predicted  vehicle 
responses  for  the  three  obstacle  tests.  In  each  case,  the  dashed 
lines  represent  the  computed  values  and  the  solid  lines  represent 
the  measured  values.  Beginning  at  the  bottom  of  each  page,  the 
first  curve  plotted  represents  simulated  and  measured  vertical 
acceleration  of  a  point  or.  the  front  axle  near  the  right  wheel, 
where  an  accelerometer  was  located  during  the  tests.  The  next 
plot  represents  vertical  acceleration  at  the  center  of  gravity  of 
the  truck  body,  and  the  third  depicts  acceleration  at  a  point  on 
the  truck  body  located  beneath  the  driver's  seat.  The  fourth 
plot  represents  computed  versus  measured  pitch  of  the  truck  body. 
It  should  be  noted  here  that  for  TesJ  No.  157  (Figure  8),  the 
pitch-measuring  instrumentation  apparently  malfunctioned  during 
the  test,  resulting  in  no  valid  record  of  measured  pitch  for  that 
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test.  At  the  top  of  each  figure,  the  position  of  the  obstacle 
relative  to  the  front  axle  of  the  vehicle  is  shown,  i.e.,  the 


obstacle  as  shown  represents  the  displacement  function  for  the 
front  tires  plotted  on  the  same  time  scale  as  the  vehicle 
responses . 

In  displaying  measured  versus  predicted  data  for  the  cross¬ 
country  run,  the  duration  of  the  run,  the  uncertainty  of  the 
exact  speed  of  the  vehicle  at  any  given  point,  and  the  possi¬ 
bility  of  yielding  or  inaccurately  defined  terrain  at  points  make 
it  impractical  to  display  actual  responses  such  as  those  for  the 
obstacle  tests.  Also  the  root  mean  square  (RMS)  of  vertical  ac¬ 
celerations  felt  by  the  driver  throughout  the  duration  of  the  test 
is  considered  to  be  of  greater  significance  than  isolated  accelera¬ 
tion  peaks.  Thus ,  the  cumulative  RMS  of  predicted  versus  measured 
vertical  acceleration  under  the  driver's  seat  for  300  feet  of  the 
cross-country  test  is  plotted  as  Figure  11.  As  with  the  previous 
figures,  the  dashed  line  represents  the  simulated  response. 
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CHAPTER  V:  DISCUSSION 
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In  comparing  the  predicted  and  measured  accelerations  for  the 
obstacle  tests (Figures  8,  9»  and  10 ),  it  should  first  be  noted  that 
the  high-frequency  oscillations  present  in  the  measured  data  are 
due  in  part  to  high-frequency  motions  such  as  frame  flexure  and 
engine  vibration  which  are  not  included  in  the  mathematical  model. 

Looking  first  at  acceleration  of  the  front  axle,  the  positive 
acceleration  produced  by  the  intial  impact  with  the  obstacle  and  the 
negative  acceleration  from  the  rebound  that  follows,  i.e.  the  first 
complete  cycle  of  axle  motion,  compare  favorably  in  shape  and  amplitude. 
This  is  especially  true  of  the  higher  speed  tests  in  which  the  ampli¬ 
tudes  of  the  acceleration  in  the  first  cycle  are  greater.  Beyond  the 
first  cycle,  the  measured  and  predicted  acceleration  curves  are  dis¬ 
placed  in  phase,  and  there  is  much  greater  measured  motion  than  simu¬ 
lated  motion.  The  axle  motion  produced  by  the  model  seems  to  have  been 
damped  out  whereas  the  measured  data  continue  to  show  oscillations. 
Generally,  the  predicted  axle  accelerations  which  represent  the  most 
sensitive  vehicle  response  that  can  be  monitored  are  better  than 
expected. 

The  center  of  gravity  (CG)  and  under  driver's  seat  (DS)  accelera¬ 
tion  curves  are  very  similar  in  all  cases — the  driver's  seat  is  located 
very  close  to  the  center  of  gravity — and  will  be  discussed  together. 

The  general  shape  of  the  predicted  curves  is  good,  particularly  during 
the  first  cycle,  but  there  is  obviously  a  problem  with  amplitude.  The 


predicted  amplitude  is  generally  too  low  throughout  the  tests, 
particularly  for  the  positive  accelerations  on  initial  impact.  The 
amplitudes  of  the  negative  accelerations  on  rebound  from  obstacle 
impact  are  much  closer  to  the  measured  data,  especially  for  Tests 
No.  157  and  162  (Figures  8  and  9»  respectively).  As  with  the  axle 
accelerations,  there  seems  to  be  some  phase  difference  beyond  the 
first  cycle.  Table  1  shows  a  comparison  of  the  simulated  versus 
measured  driver's  seat  RMS  accelerations  at  the  end  of  the  three  tests. 
While  there  is  not  close  agreement  in  the  final  RMS  values,  the  differ¬ 
ences  which  average  29  percent  for  t.he  three  tests  are  not  as  drastic  us 
the  raw  acceleration  curves  might  imply. 

In  evaluating  the  pitch  curves  (ignoring  Figure  8),  it  is  diffi¬ 
cult  to  defend  the  simulated  pitch  motion.  The  overall  trend  of 
the  pitch  curve  for  Test  No.  169  shows  some  resemblance  to  the  measured 
values  while  the  curves  for  Test  No.  l60  bear  little  resemblance. 

However,  there  must  be  some  question  as  to  the  validity  of  the  measured 
pitch  motion.  The  measured  pitch  curves  for  Tests  162  and  169  should 
be  very  similar  in  shape  and  amplitude  with  a  large  initial  negative 
pitch  and  a  positive  pitch  toward  the  end  of  the  tests  as  the  rear 
wheels  strike  the  obstacle.  There  are  considerable  differences  in 
the  two  measured  pitch  curves,  enough  to  cast  doubt  on  one  or  both 
curves.  Still,  it  must  be  admitted  that  the  model  hus  apparently 
severely  under-predicted  the  pitch  motion. 

Figure  11  shows  the  measured  versus  simulated  RMS  acceleration 
under  the  driver's  seat  for  the  cross-country  run.  The  test  was  run 
at  approximately  eight  miles  per  hour  over  rough  terrain.  The  degree 
of  agreement  between  the  measured  and  simulated  curve  is  considered  to 
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be  within  acceptable  limits.  The  average  difference  between  the 
two  curves  is  approximately  15  percent.  While  this  is  encouraging, 
the  model  is  predicting  less  than  measured  values  and  the  model 
should  ideally  predict  greater  accelerations  for  this  type  of  test. 

The  simulation  deals  with  a  nonyielding  terrain,  whereas  there  was 
probably  at  least  some  terrain  deformation  which  reduced  the  measured 
vehicle  accelerations  from  those  that  would  be  encountered  tra¬ 
versing  a  completely  rigid  riding  surface.  Also  the  straight-line 
representation  of  the  terrain  used  in  the  simulation  produces 
abrupt  changes  in  tire  deformations  that  may  be  somewhat  harsher 
than  those  produced  by  actual  terrain. 

In  trying  to  determine  reasons  for  the  lack  of  agreement  between 
the  measured  and  simulated  vehicle  responses,  particularly  the  obstacle 
tests,  there  axe  basically  two  sources  of  error  that  could  result  in 
bad  comparisons  for  an  essentially  valid  model. 

First,  the  field  test  results  are  not  repeatable.  The  forward 
velocity  of  the  vehicle  was  controlled  by  the  driver  using  the  accelera¬ 
tor  pedal.  While  he  attempted  to  maintain  constant  speeds  throughout 
the  tests,  the  speed  actually  varied  considerably  as  the  truck  ran  over 
the  obstacle.  The  simulation  maintains  a  constant  speed  absolutely. 

This-  somewhat  explains  the  apparent  phase  difference  between  the 
measured  and  computed  accelerations  beypnd  the  first  cycle.  Also  the 
steering  is  controlled  by  the  driver;  thus,  the  angle  at  which  the 
vehicle  strikes  the  obstacle  is  subject  to  human  error.  Duplicate 
field  tests  run  over  the  same  obstacle  at  the  same  speed  produce 
measured  acceleration  maximums  that  are  commonly  different  by 
fifteen  to  twenty  percent. 
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Secondly,  with  the  exception  of  the  tire  data,  the  vehicle 
parameters  used  in  the  simulation  are  questionable.  The  suspension 
parameter,  provided  by  FMC  were  measured  from  components  disassembled 
from  a  specific  vehicle,  not  the  same  vehicle  used  in  the  WES  field 
tests.  MRB  has  learned  that  characteristics  of  springs  and  shock 
absorbers  can  vary  significantly  from  unit  to  unit  of  the  same  tyre 
vehicle.  Also,  the  vehicle  moments  of  inertia  were  apparently 
computed  by  FMC  for  a  vehicle  loaded  with  a  considerably  different 
weight  distribution  than  the  vehicle  used  in  the  WES  tests.  For 
the  WES  obstacle  tests,  the  load  was  arranged  to  produce  equal  weight 
on  all  four  wheels.  Judging  from  the  pitch  motion  curves,  the  pitch 
moment  of  inertia  used  in  the  simulations  was  much  higher  than  the 
actual  moment  of  inertia  for  the  vehicle  as  loaded  for  the  obstacle 
tests . 

While  the  figures  presented  represent  the  extent  of  the  valida¬ 
tion  effort  possible  in  the  time  frame  available,  it  is  felt  that 
much  closer  comparisons  between  predicted  and  measured  data  could 
be  obtained  by  altering  the  vehicle  parameters,  i.e.  "tuning"  the 
model  to  match  the  measured  accelerations.  However,  this  would  be  a 
time-consuming  process  and  would  not  necessarily  be  a  valid  approach 
to  the  problem.  Ideally,  complete  validation  should  be  accomplished 
using  vehicle  parameters  obtained  through  tests  and  measurements  on 
the  field  test  vehicle,  assuring  parameter  reliability. 
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In  an  attempt  to  establish  criteria  for  control  of  the  integration 
step  size  for  the  model,  simulation  of  obstacle  traversal  was  executed 
using  several  different  constant  time  intervals.  First,  using  a  forward 
velocity  corresponding  to  the  highest  speed  validation  test  (Mo.  169), 
a  "reference"  time-step  was  established  by  executing  the  complete 
simulation  at  a  constant  time-step,  reducing  the  time-step  by  half  and 
executing  again.  This  was  repeated  until  further  reduction  in  the  time- 
step  produced  no  detectable  change  in  the  axle  acceleration  curve,  axle 
acceleration  representing  the  most  sensitive  vehicle  response.  Ex¬ 
ecution  of  the  simulation  at  this  reference  time-step  represents  *  he 
"best"  simulation  results  possible  from  the  model. 

Having  established  the  reference  time-step  (0.00199  sec.  in  this 
case),  plots  were  made  to  compare  vehicle  responses  produced  by  the 
reference  simulation  with  responses  produced  by  simulations  executed 
at  time-steps  of  (2n)  T,  where  T  is  the  reference  time-step  and  n  is  a 
positive  integer.  Figures  12  through  15  show  the  results  of  this  pro¬ 
cedure.  The  responses  chosen  for  comparison  were  front-axle  accelera¬ 
tion,  velocity,  and  displacement,  and  acceleration  under  the  driver's 
seat.  Axle  acceleration  represents  the  response  in  which  effects  of 
changes  in  time-step  are  most  obvious;  axle  velocity  and  displacement 
are,  in  effect,  the  forcing  functions  to  the  vehicle  body;  and  the 
driver's  seat  acceleration  is  actually  the  vehicle  response  of  primary 

concern.  In  each  plot,  the  solid  lines  reflect  the  reference  time- 
step  responses  and  the  dashed  lines  the  simulation  at  the  larger  time- 
step.  Also  plotted  with  the  axle  velocity  and  displacement  are  the 


is. oi  nos iimusjs  "'0  mu  ;k3vi3’'n„<:  :  ji«a  ois/mii  i  men  3  a  jiiu 


BCCELERRTJON  (G*Sl 


15.0!  N0I14«91333« 


>■ 


- 2. _ 


error  terms,  computed  by  Merson's  error  formula,  for  these  two 
integrations.  These  are  plotted  as  discrete  points,  the  scales  for 
these  points  being  drawn  on  the  right  side  of  the  figure. 

Beginning  with  Figure  12,  the  dashed  lines  represent  simulation 
at  double  the  reference  time-step.  While  the  axle  acceleration 
shows  some  departure  from  the  reference,  all  the  other  curves  are 
essentially  the  same.  Obviously,  this  time-step  would  be  acceptable 
throughout  the  simulation.  Figure  13  reflects  use  of  a  time-step 
of  four  times  the  reference.  Here  the  axle  acceleration  shows  a 
further  departure  from  the  reference,  and  the  velocity  and  displac:-- 
rcent  curves  are  beginning  to  show  a  slight  deviation.  However,  the 
driver's  seat  acceleration  is  still  essentially  unchanged.  Looking 
at  Figure  l1*,  which  represents  a  time-step  of  eight  times  the 
reference,  the  uxle  acceleration  has  departed  signi f icantly  from  the 
reference  curve,  the  velocity  and  displacement  show  some  error, 
but  the  driver's  seat  acceleration  curve  is  still  very  close  to 
the  reference  curve.  Doubling  the  time-step  again  (now  sixteen 
tires  the  reference)  yields  Figure  j.5.  Here  the  axle  acceleration 
is  severely  degraded  and  the  velocity,  displacement,  and  driver's 
seat  acceleration  curves  show  u  further  departure  from  the  reference 
than  did  curves  produced  by  smaller  time-steps,  but  the  critical 
driver's  scat  acceleration  still  shews  very  good  shape  and  amplitude 
and  is  still,  for  all  practical  purposes,  acceptable  for  computing 
RMS  acceleration.  Attempting  to  increase  the  step-size  to  32T 
produced  a  condition  in  which  the  maximum  value  of  one  of  the 
suspension  spring  force-deflection  tables  was  exceeded  (a  condition 
described  in  Chapter  III  as  indicating  need  for  time-step  reduction). 
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thus  making  further  simulations  at  larger  constant  time  intervals 
inpracticable  for  this  vehicle  speed. 

Hot  being  able  to  Judge  as  unacceptable  the  results  from  any 
of  the  time  steps  used  for  the  22  ft. /sec.  simulation,  it  was 
decided  to  repeat  the  procedure  above,  simulating  obstacle  traversal 
at  twice  the  forward  velocity,  i.e.  Mi  ft. /sec.  'flic  reference 
time-step  was  established  ns  0.30093  seconds  using  the  same  procedure 
as  previously  described.  Figure  l6  shows  the  results  of  simulation 
at  twice  the  reference,  obviously  revealing  no  significant  deviation 
in  any  of  the  responses.  Figures  IT,  18,  and  19  reflect  the 
responses  computed  using  time-steps  of  hT ,  8T,  and  l6T,  respectively. 
These  reveal  essentially  the  same  pattern  as  would  be  expected  based 
on  the  previous  series  of  figures,  i.e.,  some  departure  from  the 
reference  axle  accelerations,  velocities,  and  displacement,  but 
little  significant  difference  in  the  driver's  seat  acceleration 
curves.  However,  Figure  20,  which  represents  execution  using  an 
integration  time-step  of  32T,  deserves  further  comment.  Here  the 
axle  acceleration  is  severely  degraded  bearing  little  resemblance  to 
the  reference  curve.  With  this  type  of  axle  acceleration  curve,  one 
would  expect  that  the  vehicle  responses  in  general,  including 
driver's  seat  acceleration,  would  show  a  significant  departure 
from  the  reference  curves.  This  is  not  the  case  here.  Axle 
velocity  and  displacement  and  consequently  the  driver's  seat  accel¬ 
eration  still  show  the  essential  effects  of  traversal  of  the 
obstacle.  For  the  driver's  seat  acceleration,  the  positive 
amplitude,  the  time  of  the  positive  peak,  the  negative  amplitude, 


and  the  time  of  the  negative  peak  are  all  close  to  the  reference 
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values.  Attempting  to  execute  the  simulation  at  a  time-step  of 
6  b  T  produced  the  condition  described  earlier  in  which  a  vehicle 
parameter  table  was  exceeded.  Again,  it  was  not  possible,  within 
the  realistic  limits  of  suspension  components,  to  execute  simulation 
at.  a  time-step  large  enough  to  severely  degrade  the  vertical  accel¬ 
eration  at  the  significant  point  on  the  vehicle  body,  i.e.  the 
driver's  seat. 

To  aid  in  extracting  conclusions  from  the  figures,  Table  2 
represents  a  summary  of  the  results.  Column  1  specifies  the  forward 
velocity  in  ft. /sec.  (22  or  UU ) .  Column  2  is  the  time-step  being 
compared  with  the  reference.  Column  3,  labeled  "AS,"  indicates  the 
forward  motion  (inches) of  the  vehicle  per  time-step.  Column  b  is 
the  maximum  error  (in. /sec.)  computed  by  the  Merson  formula  in  the 
integration  yielding  axle  velocity.  Column  5  is  the  maximum  error 
(inches)  computed  by  the  Merson  formula  in  the  integration  yielding 
axle  displacement.  Column  6,  labeled  "RMS  Acc.  Residual,"  is  the  abso¬ 
lute  difference  (G's)  between  the  compared  and  the  reference  driver's 
seat  RMS  acceleration  at  the  end  of  the  simulation.  Column  7> 
labeled  "Avg.  Acc.  Residual,"  represents  the  average  absolute 
difference  (G's)  between  the  compared  and  the  reference  driver's 
seat  accelerations  for  all  discrete  points  in  the  simulation. 

How  can  these  results  be  related  to  the  basic  questions  con¬ 
cerning  time-step  management  posed  in  Chapter  III?  First,  consider 
determination  of  maximum  time-step.  As  stated  previously,  it  is  not 
uncommon  to  relate  maximum  time-step  to  the  period  of  vibration  of 
the  highest  frequency  component  being  modeled.  In  this  case  this 
means  vertical  motion  of  the  axle,  which  for  the  vehicle  simulated 
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vibrates  at  approximately  ten  cycles  per  second.  This  would 
indicate  a  maximum  time-step,  using  ten  integrations  per  period,  of 
about  0.01  second.  However,  in  simulations  at  both  22  ft. /sec. 
and  hi*  ft. /sec.,  time-steps  in  excess  of  0.01  second  produced  ac¬ 
ceptable  results.  In  fact,  it  was  seen  that  the  accuracy  of  the 
modeled  axle  acceleration  was  not  critical  to  producing  acceptable 
vehicle  body  accelerations.  To  confirm  the  apparent  conclusion  that 
the  time-step  need  not  be  limited  by  1/10  the  axle  vibration  period, 
the  cross-country  run  described  previously  was  simulated  at  a 
time-step  of  0.03125  second  which  is  approximately  three  times  this 
limitation.  The  driver's  seat  RMS  acceleration  history  produced  is 
virtually  identical  with  that  plotted  on  Figure  11,  which  was  produced 
using  a  time-step  of  0.00781  second.  This  indicates  that  the  period 
of  axle  vibration  warrants  consideration  when  choosing  maximum 
time-step  only  if  the  axle  acceleration  itself  is  of  primary  impor¬ 
tance,  which  is  not  the  case  for  the  current  study.  Thus,  the 
maximum  time-step  can,  in  general,  be  limited  by  a  maximum  allowable 
forward  displacement  of  the  vehicle  in  one  time-step.  To  avoid 
missing  significant  terrain  features  or  creating  unrealistic  tire 
deflections  by  too  great  a  forward  advance,  the  time-step  should 
be  limited  such  that  the  forward  displacement  of  the  vehicle  in 
one  t-L.-e-step  does  not  exceed  the  diameter  of  the  tire  body.  Note 
that  the  last  step  size  used  (Mi  ft. /sec.)  exceeds  this  limitation 
and  would  not  be  used  in  the  production  computer  program  (in  which 
this  maximum  interval  criterion  is  now  implemented.) 

As  stated  in  Chapter  III,  the  minimum  allowable  time-step  can 
be  established  to  correspond  to  some  minimum  forward  displacement, 
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i.  e.,  that  forward  displacement  corresponding  to  the  "reference" 
time-step  defined  previously.  Both  reference  time-steps  established 
(0.00195  at  22  ft. /sec.  and  0.00098  at  Ult  ft. /sec.)  correspond  to 
forward  displacement  of  approximately  0.5  inch.  It  seems  reasonable 
then  to  assume  that  the  minimum  time-step  can  be  established  as  a 
function  of  forward  velocity,  i.e.,  that  time-step  that  corresponds 
to  forward  displacement  of  0.5  inch.  Attempts  to  reduce  the  time- 
step  below  this  value  will  produce  termination  of  the  simulation. 

Turning  next  to  the  questions  regarding  time-step  reductions, 
if  driver's  seat  acceleration  is  the  only  time  history  of  Interest, 
then  the  computed  Merson  error  is  not  as  significant  a  factor  in 
time-step  reduction  as  previously  thought.  Operating  at  the  maxi¬ 
mum  time-step  as  specified  above,  reducing  the  time-step  only  when 
the  requirement  is  indicated  by  excessive  deflections  or  velocities 
or  by  geometric  considerations  would  probably  produce  acceptable 
results  in  most  cases.  However,  it  is  wise  to  monitor  the  Merson 
error,  as  the  other  conditions  mentioned  are  not  completely  reliable 
indicators  of  time-step  problems,  e.g.,  a  forward  step  too  large  for 
current  terrain  conditions  does  not  always  produce  a  table-look-up 
error.  In  determining  which  integration  error  to  monitor,  it  was 
recognized  that  if  axle  acceleration  is  tne  vehicle  response  most 
sensitive  to  changes  in  terrain  conditions,  it  follows  that  inte¬ 
gration  of  axle  accelerations  to  produce  velocities  should  provide 
the  raosn  responsive  error  indicator.  It  was  decided  that  velocity 
errors  should  be  evaluated  on  an  absolute  rather  than  a  relative 
basis  since  relative  error  evaluation  produces  the  potential  of 
over-rating  the  error  significance  at  points  where  the  magnitude 


of  the  integration  variable  (velocity)  is  small,  e.g.,  near  zero 
crossing.  From  Table  2  it  appears  that  the  computed  error  does 
somewhat  reflect  degradation  in  vehicle  responses.  However,  the 
computed  error  is  obviously  not  a  function  of  h^;  doubling  the 
time-step  has  not  produced  error  terms  increased  by  a  factor  of 
32.  In  one  case,  doubling  the  time-step  actually  resulted  in 
decreased  maximum  reported  error.  Moreover,  the  computed  error  does 
not  correlate  well  with  the  RMS  acceleration  residual  or  the  average 
acceleration  residual.  Thus,  the  determination  of  a  maximum  allow¬ 
able  integration  error  must  be  made  on  a  rather  arbitrary  basis. 

From  Figures  12  through  20  and  the  corresponding  maximum  velocity 
errors,  it  seems  acceptable  to  establish  a  maximum  allowable  inte¬ 
gration  error  at  10  in. /sec.  This,  in  conjunction  with  the  other 
conditions  that  produce  step  size  reduction,  should  effectively  reduc 
the  time-step  when  required.  The  role  of  time-step  reduction  assumes 
much  less  significance  when  a  realistic  maximum  time-step  has  been 
established  as  previously  defined. 

Finally,  considering  possibilities  for  time-step  increase,  it  is 
apparent  from  the  figures  that  curves  that  very  closely  approximate 
the  reference  are  accompanied  by  very  small  integration  errors.  Thus 
increase  in  time-step  based  on  computed  integration  error  is  appar¬ 
ently  a  valid  procedure  for  this  model.  Again,  as  with  the  maximum 
error,  the  minimum  error  selection  is  an  arbitrary  decision  since 
there  is  no  well-defined  relationship  between  integration  error  and 
degradation  of  vehicle  response.  However,  an  arbitrary  choice  will 
accomplish  the  objective  which  is  to  increase  the  step-size  back 
to  the  maximum  as  soon  as  possible  after  a  severe  disturbance 


(such  as  the  obstacle  used  here)  has  been  passed.  Based  on  Figures 

12  through  20,  an  acceptable  minimum  error  can  be  established  at 

1  in. /sec.  A  computed  axle  velocity  error  less  than  1  in. /sec.  will 

result  in  the  step  size  being  doubled  for  the  next  integration  step. 

\ 

It  is  interesting  to  note  the  small  peaks  present  on  the 
axle  acceleration  curves,  produced  using  small  time  intervals,  par¬ 
ticularly  those  produced  by  the  reference  curves.  These  small 
abrupt  changes  in  axle  acceleration  can  be  attributed  to  discontin¬ 
uities  in  tire  stiffness  inherent  in  the  segmentation  scheme 
described  previously.  Each  small  peak  can  be  related  to  a  segment 
making  or  breaking  contact  with  the  obstacle  as  the  tire  moves 
forward.  While  undesirable  effects  of  segmentation  are  apparent 
when  viewing  the  axle  acceleration  curves,  the  integration  process 
tends  to  have  a  smoothing  effect  (note  the  velocity  and  displacement 
curves),  resulting  in  no  adverse  effeci  s  being  "felt"  by  the  vehicle 
body. 


Computing  Cost 

Cost  of  using  the  model  is  directly  proportional  to  the  number 
of  time-steps  required  for  the  simulation.  On  the  Honeywell  Ghho 
the  computing  time  required  is  approximately  two  seconds  central 
processor  (CP)  time  per  time-step.  This  means,  for  simulation  of 
typical  cross-country  traversal,  about  300  to  500  seconds  computing 
time  per  100  feet  of  terrain.  A  typical  charge  for  CP  time  for  this 
class  of  computer  is  about  $.03  per  second.  Thus  cost  for  simulation 
of  a  1000-foot  cross -country  would  be  approximately  $90  to  $150. 
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A  "number  crunching"  program  such  as  this  would  normally  be  much 
more  economical  to  run  on  larger,  faster  computing  systems. 
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CHAPTER  VI:  CONCLUSIONS  AND  RECOMMENDATIONS 

Conclusions 


The  results  of  this  study  suggest  the  following  conclusions: 

A.  Vehicle  vibrations  can  be  predicted  within  acceptable  error 
limits  using  a  numerical  model  (digital  simulation)  that  describes 
the  vehicle  as  a  rigid  body  connected  to  rigid  axles  by  springs  and 
viscous  dampers. 

B.  The  presence  of  discontinuities  in  the  representation  of 
riding  surface  profiles,  tire  flexure,  and  suspension  springs  and 
dampers  does  not  significantly  degrade  the  performance  of  the  model. 

C.  The  Runge-Kutta-Merson  numerical  integrat.' on  method  is  an 
effective  technique  for  integrating  the  differential  equations 
describing  the  vibratory  motions. 

D.  The  numerical  integration  interval  can  be  dynamically  con¬ 
trolled,  based  on  simulation  responses  and  computed  integration 
errors,  to  provide  improved  model  performance. 

1.  The  period  of  axle  vibration  should  be  considered  when 

establishing  a  maximum  time-step  only  if  axle  accel¬ 
eration  itself  is  of  primary  importance. 

2.  Determination  of  vehic]  iy  accelerations  requires 

that  the  maximum  time-step  be  limited  such  that  the 
forward  displacement  of  the  vehicle  in  one  time-step 
does  not  exceed  the  tire  cross-section  height. 

3.  The  minimum  time-step  is  that  time-step  that  corresponds 

to  forward  displacement  of  0.5  inch. 

i».  Simulation  responses  that  yield  displacements  or  veloc¬ 
ities  in  excess  of  physical  limits  of  the  system 
require  a  time-step  reduction. 

5 .  Integration  error  computed  by  the  Merson  formula  can 

be  used  for  the’  purpose  of  time-step  reduction.  An 
error  of  10  in. sec.  in  vertical  velocity  of  any  axle 
can  be  used  as  a  maximum  allowable  integration  error. 

6.  Integration  errors  computed  by  the  Merson  formula  can  be 

used  for  the  purpose  of  increasing  the  time-step. 
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If  velocity  errors  for  all  axles  are  less  than 
1  in. /sec.,  a  time-step  increase  is  justified. 

7.  The  initial  time-step  she  Id  be  the  same  as  the  maxi¬ 
mum  time-step. 


Recommendations 


It  is  recommended  that 

A.  The  validation  of  the  model  be  completed  by 

1.  Precise  determination  of  vehicle  parameters  for  a  test 

vehicle . 

2.  Execution  of  field  tests  so  that  the  variables  of  the 

tests,  i.e.  forward  velocity,  angle  of  approach,  etc., 
are  more  precisely  controlled,  thus  eliminating  the  driver 
response  from  the  tests. 

3.  Simulation  of  the  field  tests  and  comparison  of  the 

measured  and  simulated  vehicle  responses. 

B.  A  study  be  made  to  determine  the  minimum  number  of  tire 
segments  that  represents  tire  flexure  to  the  degree  of  accuracy 
necessary  to  yield  acceptable  vehicle  body  accelerations. 

C.  The  model  be  extended  to  include  simulation  of  vehicles 
having  other  types  of  suspension  systems,  e.g.,  independent  sus¬ 
pensions  and  "walking  beam"  suspensions. 


ABSTRACT 


Windell  Francis  Ingram,  Master  of  Science,  1972 
Major:  Engineering  Mechanics,  College  of  Engineering 
Title  of  Thesis:  A  Numerical  Model  of  the  Ride  Dynamics  of  a 
Vehicle  Using  a  Segmented  Tire  Concept 
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ABSTRACT 


The  purpose  of  this  study  vas  to  develop  and  validate  a 
digital  computer  simulation  of  the  ride  dynamics,  including  bounce, 
pitch,  and  roll,  of  a  general  single-hull,  solid-axle  vehicle  with 
an  arbitrary  number  of  axles.  The  study  also  included  an  effort 
to  minimize  the  execution  time  and  thus  the  cost  of  the  simulation 
by  providing  automatic  determination  of  the  optimum  numerical  inte¬ 
gration  interval  to  be  used  at  each  point  in  the  simulation. 

The  vehicle  was  defined  as  a  rigid  body  connected  to  rigid  axles 
by  springs  and  viscous  dampers.  Tires  were  described  as  clusters 
of  radially  projecting  springs  being  deflected  by  a  nonyielding 
riding  surface.  The  suspension  springs  and  dampers  were  described 
by  force-deflection  and  force-rate  tables,  thus  allowing  nonlinear 
characteristics  to  be  handled  by  a  piece-wise  linear  representation 
with  linear  interpolation  between  the  points. 

The  numerical  integration  method  used  was  the  Runge-Kutta-Merson 
method,  which  is  a  fourth-order  method.  The  method  requires  five 
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derivative  evaluations  per  integration  and  provides  an  estimate  of 
the  truncation  error  that  can  he  used  for  automatic  interval 
adjustment. 

Field  tests  of  a  representative  vehicle  traversing  a  single 
obstacle  at  three  different  speeas  and  a  cross-country  terrain  at 
one  speed  were  simulated.  Graphs  are  presented  comparing  simulated 
vehicle  responses  with  field  test  data  provided  by  the  Waterways 
Experiment  Station. 

The  effects  of  integration  step-size  on  simulated  vehicle 
responses  are  illustrated  with  graphs  comparing  simulated  vehicle 
responses  using  several  different  integration  intervals.  A  dis¬ 
cussion  of  criteria  for  automatic  management  of  the  integration 
interval  is  presented. 
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generalized  vehicle  dynamics  PROGRAM  for  single 
HULLiSOLID  AXLE  VEHICLES  ;  USES  A  SEGMENTED  TIRE 
CONCEPT  ,  THIS  PROGRAM  Is  A  MODIFICATION  OF  A 
PROGRAM  DEVELOPED  by  fmc  for  the  waterways  experiment 
station. 


difinitions  of  variables 


NAX  NUMBER  OF  AXLES 

vel  forward  velocity 

accel  forward  acceleration 

BETA  TRAVEL  PATH  ANGLE 

alpha  terrain  slope  angle 

deltim  integration  step  size  for  equilibrium  computation 

TSTOP  STOP  TIME 

HMIN  minimum  STEP  SIZE  FOR  INTEGRATION 

HSTOP  MAXIMUM  TRAVEL  DISTANCE  FOR  VEHICLE 

ACGHT  HEIGHT  OF  AXLE  CG  AT  START  OF  EQUILIBRIUM  COMPUTATIO 

BCGHt  HEIGHT  OF  BODy  CG  AT  START  OF  EQUILIBRIUM  COMPUTATIO 

LAX(I)  DISTANCE  FROM  BODY  CG  TO  AXLE 

LSUS(J)  DISTANCE  FROM  BODY  CG  TO  SUSPENSION 

LASUS(J)  DISTANCE  FROM  AXLE  CG  TQ  SUSPENSION 

LTIR(J)  DISTANCE  FROM  AXLE  CG  TO  TIRE 

WTB  BODY  HEIGHT 

WTAX(I)  AXLE  WEIGHTS 

IBP  BODY  PITCH  MOMENT  OF  INERTIA 

ibr  body  roll  moment  of  inertia 

IAXR  axle  ROLL  MOMENTS  OF  INERTIA 

KSAME(J)  SAMENESS  TABLE  FOR  SUSPENSIONS  AND  TlRES 

KSTN  SUSPENSION-TIRE  NUMBER 

KSMN  SAMENESS  NUMBER 

WRAD  WHEEL  RADIUS 

SREFD(J)  SPRING  REFERENCE  DISTANCES 

TREFO(J)  TIRE  REFERENCE  DISTANCES 

T ( K  )  TEMPORARY  STORAGE  FOR  CARD  DATA 

NSEG  NUMBER  OF  TIRE  SEGMENTS  (MAXIMUM  24) 

DSEG(I)  SEGMENT  CENTERLINE  W.R.T,  VERTICAL  (DEGREES) 

SEGK(l)  SEGMENT  SPRING  CONSTANT 

deltaik, j> ,fs<k, j)  profiles  of  suspension  spring  force 

DELTD(K,J) ,FD(K, J)  PROFILES  OF  SUSPENSION  DAMPING  FORCE 

DELDtK. J) ,FFD(K, J)  PROFILES  OF  TIRE  DAMPING  FORCE 

WTOT  TOTAL  WEIGHT  OF  TRUCK 

FDR  I  V  HORIZONTAL  PROPULSION  FORCE 

WLAG(I)  LAG  DISTANCES  OF  each  WHEEL  CENTER  BEHIND  LEAD  W 

HH  INTEGRATION  STEP  SIZE 

DU)  XDDB  ACCELERATION  OF  BODY  CG 

D  (  2  )  XDB  VELOCITY  OF  BODY  CG 

D ( 3 )  THDDB  BODY  PITCH  ANGLE  ACCELERATION 

0(4)  THDB  BODY  PITCH  ANGLE  VELOCITY 
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•  D(5)  GDDB  BODY  ROLL  ANGLE  ACCELERATION 

•  D ( 6 )  GOB  BODY  ROLL  ANGLE  VELOCITY 

•  D  ( 7  )  D  ( 16  )  XDDA(I)  ACCELERATION  OF  AXLE  CG 

•  D < 1 7 ) -D { 26 )  X D A ( I )  VELOCITY  OF  AXLE  CG 

•  D ( 2 7 ) -D ( 36 )  PHDDA(I)  AXLE  ROLL  ANGLE  ACCELERATION 

•  D ( 37 ) -D ( 46  )  PHD A ( 1 )  AXLE  ROLL  ANGLE  VELOCITY 

•  D ( 4 7 )  HDD  HORIZONTAL  ACCELERATION  OF  TRUCK 

•  D <  46 )  HD  HORIZONTAL  VELOCITY  OF  TRUCK 

•  D < 49 )  WORKD  TIKE  DERIVATIVE  OF  WORK 

•  D ( 50 )  XLDD  COHPUTED  HORIZONTAL  VELOCITY  AT  CG 

o  D ( 51 )  TDOT  TIME  DERIVATIVE  OF  TIME 


V  ( 1 ) 

XDB 

VELOCITY  OF  BODY  CG 

V  ( 2  ) 

XB 

POSITION  OF  BODY  CG 

V  ( 3 ) 

thdb 

BODY  PITCH  ANGLE  VELOCITY 

V  (  4  ) 

THB 

BODY  PITCH  ANGLE 

V  ( 5 ) 

GDB 

BODY  ROLL  ANGLE  VELOCITY 

V  (6) 

GB 

BODY  ROLL  ANGLE 

V ( 7 ) -V ( 16  ) 

XDA  ( I  ) 

VELOCITY  OF  AXLE  CG 

V l l7 ) -V { 26 ) 

XA  (  I  ) 

POSITION  OF  AXLE  CG 

V ( 27 ) -V  <  36  ) 

PHDA  ( I ) 

AXLE  ROLL  ANGLE  VELOCITY 

V ( 37  ) - V ( 46  > 

PHA  <  I ) 

AXLE  ROLL  ANGLE 

V  (  4  7  ) 

HD 

HORIZONTAL  VELOCITY  OF  TRUCK 

V  (  46  ) 

HP 

HORIZONTAL  POSITION  OF  LEAD  WHEEL  CENTE 

V  (  4  9  ) 

WORK 

WORK  (INCH-POUNDS) 

V(50> 

SB 

NOT  USED 

V  { 51 ) 

TIME 

CURRENT  TIME  (SECONDS) 

NSUS  NUMBER  OF  SUSPENSIONS 

DLS(J)  DEFLECTIONS  OF  SUSPENSIONS 

DLSD(J)  DEFLECTION  VELOCITIES  OF  SUSPENSIONS 

DLTDCJ)  DEFLECTION  VELOCITIES  OF  TIRES 

SPS(J)  SPRING  FORCES  OF  SUSPENSIONS 

DMS(J>  DAMPING  FORCES  OF  SUSPENSIONS 

DMT(J)  DAMPING  FORCES  OF  TIRES 

►*  <  J  >  TOTAL  SUSPENSION  FORCES 

FF ; J )  TOTAL  TIRE  FORCES 

SUM0  SUM  OF  FORCES  (TEMPORARY) 

SUMOM  SUM  OF  MOMENTS  (TEMPORARY) 

FHORIZ  SUM  OF  HORIZONTAL  FORCES 

HORMOM  SUM  OF  HORIZONTAL  MOMENTS 


FFV(I)  TIRE  FORCE  VECTOR  SUM 

XWCN  HORIZONTAL  POSITION  OF  WHEEL  CENTER 

XRET  REAR  EXTREME  OF  TIRE  (HORIZONTAL) 

XF  t  T  FRONT  EXTREME  OF  TIRE  (HORIZONTAL) 

CC(l)  POSITION  VECTOR  OF  WHEEL  CENTER 

TMAG  MAGNITUDE  OF  RXY  VECTORS 

A A ( I )  TEMPORARY  STORAGE 

ST  SLOPE  OF  TERRAIN 

sm  slope  of  wheel  segment 

XI, Yl  COORDINATES  OF  HUB 

X2.Y2  COORDINATES  OF  NEAR  TERRAIN  STATION 

X3,Y3  COORDINATES  OF  FAR  TERRAIN  STATION 

RXY(I.J)  POINTS  ON  WHEEL  CIRCUMFERENCE 


CX ( I , J )  component  OF  SEGHENT  centerline 
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«  CY  <  J  ,  J  ) 

•  V  V  V  (  I  ) 

•  VELV  ( I ) 

•  sdld 

•  TSFOR 

®  FORVeC( I ) 
»  FVFR(J) 

»  flat ( j ) 


COMPONENT  OF  SEGMENT  CENTERLINE 
UNIT  VECTOR  IN  direction  OF  TIRE  FORCE 
wheel  CENTER  VELOCITY  VECTOR 
TIRE  DEFLECTION  VELOCITY 

tire  spring  force 

TIRE  FORCE  VECTOR 

VERTICAL  COMPONENT  OF  TOTAL  SUSPENSION  FORCE 
LATERAL  COMPONENT  OF  TOTAL  SUSPENSION  FORCE 
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REAL  NAX,  NPPNT,  LSUS.  LTIR,  lax,  IBP,  IBR,  IAXR,  LaSUS 

DIMENSIONS  FOR  THE  FOLLOWING  variables  allow. 

A  1U  axle  VEHICLE  .  THESE  DIMENSIONS  MUST  NOT  BE 
CHANGED  UNLESS  THE  RELATED  EQUIVALENCE  STATEMENTS 
AND  THE  RELATED  SUBSCRIPTING  IN  THE  DER I  V  I T I VE 
ROUTINE  ARE  CHANGED  ALSC. 

COMMON  0(51),  V  (  5  1  )  >  VL  <  51 ) /  P1<51)»  P2(51), 

P3 ( 51  )  )  P4(51),  P  5 ( 5 1 ) 

opOOo'>oo#r<»ooofl'>»»<»* 

, . .DIMENSIONS  For  THE  FOLLOWING  VARIABLES  ALLOW 
A  3  AXlE  VEHICLE.  THEY  WERE  REDUCED  TO  CONSERVE 
CORfc  SPACE.  THESE  DIMENSIONS  MAY  BE  CM ANGE D ( F OR 
IIP  TO  10  AXLES '/WITHOUT  MAKING  AMY  OTHER  CHANGES, 

Common  lax<3),  lsuS(6)i  l t I R < 6 ) .  wTAX(3>, 

I  A  x  r  (  3  )  ,  CELTA(?5,6),  FS ! 25 , 6 ) ,  DElTD<25,6), 

F  0 ( 2  5  »  R  )  »  DSE  G ( 25 / 6  )  ,  SEGK(25,6>;  DELD(25,6>» 

FFD ( 25  <  6  )  .  K  S  A  M  E ( 6  )  .  DLS ( 6  > ,  DLSD<6>,  DLT(6>, 

DL  T  D ( 6  )  «  SPS ( 6  )  .  DMS  (  6  )  ,  SpT<6),  DHT(6), 

F ( 6 ) ,  FF( 6),  F  vER ( 6  )  i  FLAT<6)»  SREFD(6), 

TREFD(R).  N  S  E  G  (  6  )  »  L  A  SUS ( 5 ) >  WLAG<6>,  CX(25,6>,  C Y  (  25 . 6 ) 


...DIMENSIONS  on  REMAINING  VARIABLES  ARE  NOT 

l  function  of  the  maximum  number  of  allowable 
i.  X  L  f.  S  , 

common  a  a  c  2  > ,  CC(2>.  vgl  v  <  2 ) ,  vw<2>#  fqrvec<2), 

F  F  V ( 2  )  j  RXY ( 2 i 25  )  i  TMAG!25) 

common  icomnt(24>,  fhoriz,  xwcn.  xrit,  xfet, 

SDI.D,  Y  3  r  OR.  TDFOR,  SUMFOR,  NINVARi  MOTION, 
NAX,  1 BR ,  BP,  HH.  ITERR,  ITERL,  FDRIV, 

ACCEL.  WTOT,  WTB.  VEL/  GG,  LUERR,  HMIN, 

HMAX,  NmN,  JS/  I T  R ,  WRAD,  IWHL,  USJ,  JAJ, 

1 NDx ,  FlTOUT ,  JAX ,  NSUS 

Common  s  t  a  c  5  o ) ,  hgt<50) 
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c 

DIMENSION  T  (  8  ) »  NMFL  ( 2 ) 

C 

DIMENSION  XDDA(IO),  XDA(10>,  XAUO),  PHDDAUO), 

4  PHDA(lO),  PHA(IO) 

c 

EQUIVALENCE  (XDDB.  D(l))#  (TKDDB,  DC3>>, 

4  (  GDDB  ,  D  ( 5  )  ) .  (XDDAU),  D  (  7  >  > »  (PHDDA(1>, 

4  D ( 27 ) ) t  (HDD.  D ( A  7  ) ) »  (SDDB.  D(49>),  (TDOT,  D(5l>) 

equivalence  <xdb,  van.  <xb,  v(2>).  (thdb, 

4  V ( 3  )  )  ,  ( THB i  V ( A  )  ) ,  (  GDB .  V(5>>,  (GB.  V  t  6  )  )  » 

4  (XDA(l).  V(7)),  ( X A ( 1 ) .  V  < 1 7  > ) ,  (PHDA<1), 

4  V ( 27 )  ) »  (PHA(l),  V(37)>,  ( HDi  V(47)>,  (HP, 

4  V ( 48  > ) .  (SDB,  v ( 49  ) ) i  (SB.  V(50>>.  (TIME,  V(5l)> 

C 

EQUIVALENCE  (PLTOUT,  NMFL> 

C 

loo  format  ( seio . o ) 
no  format  (1615) 

120  FORMAT  CO".  4X,  "EQUILIBRIUM  CONDITIONS"  //  (5X,  9F12 , 7 )  ) 
...COMMENTS  hill  BE  WRITTEN  TO  THE  OUTPUT  FILE 
PRINT,  "COMMENTS- 

Read  130,  icomnt 
130  format  (2aa3> 

. . ,  input  file  name  is  PARAMS 

call  openfu,  -params-) 

. . .ESTABLISH  output  file 

print,  -output  file  name-,  *  * 
read  140,  pltout 
14 u  Format  { a 6 ) 

write  (3)  dummy 
call  closef«3.  pltout) 
call  openf(3,  pltout) 

PRINT,  "SPEED  OF  TRAVERSAL  IN  IN/SEC",  ♦  • 

READ,  VEL 
Gq  =  386.064 

DtR  «  3.141592653  /  100-0 

read  input  data  ••••***• 

*»  Ov  >»•.#*»«#»*»»•«• 


motion  =  o 

...  read  constants  (Dummy  is  inserted  for  variable 
fields  no  longer  used, ) 

Read  (l;  100)  DUMMY.  ACCEL,  BETA,  ALPHA 

READ  (1;  100)  DELTIM 

READ  (1;  100)  NAx.  ACGHT,  BCGHT 
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1 . 


nno 


Read  ti;  loin  wtej,  ibp,  i br 

JAX  a  NAX 
NSUS  =  2.0  *  NAX 

C  READ  AXLE  DISTANCES,  suspension  and  tire  distances 

read  (1;  100 )  ( L AX (  I  ),  1=1,  JAX) 
read  (IS  1 0  0 )  (LSUS(l),  1  =  1.  NSUS) 

READ  (1;  100)  ( L ASUS  < I  )  >  I  =  1,  NSUS) 

Read  a;  100)  (ltiru),  i  =  1,  nsus) 

C  READ  axle  WEIGHTS  AND  MOMENTS 

Read  a;  iooj  <wtax<i),  i  =  1,  j a x ) 

read  (1;  100)  (IAXR(l),  I  =  1.  JAX) 

C  ...  READ  REFERENCE  DISTANCES 

READ  (l;  100)  (SREFD(l),  i  =  1,  NSUS) 

READ  (1,  100)  (TREFD(I),  I  =  1.  NSUS) 

WRAD  =  TREFD(l) 

C  ...READ  NUMBER  OF  SEGMENTS  PER  TIRE 

Read  <i;  no)  <nseg(I),  i  =  1,  nsus) 

...  input  force  profile  data  ... 

. . .read  parameter  card 
150  read  (i;  no)  kstn,  ksmn 

c  ...  test  for  flag 

If  (KSTN)  330.  330.  160 
C  ...  SET  SAMENESS  NUMBER 

160  KSAME(KSTN)  =  KSMN 
C  ...  TEST  FOH  PROFILE  DATA 

IF  (KSMN)  150.  170.  150 

170  continue 

C  ...  READ  SUSPENSION  SPRING  FORCE  PROFILE 

It  =  l 

iso  read  a;  100)  t 

IF  ( T ( 1  )  ♦  T ( 3  )  )  190,  210.  190 
190  Do  20C  I  =  1,  7,  2 

DELTA< IT, KSTN)  =  T(I  ) 

F  S ( IT.KSTN)  =  T  (  I  ♦  1  ) 

IT  =  IT  ♦  1 

200  continue 

Go  TO  180 

C  ...  READ  SUSPENSION  DAMPING  FORCE  PROFILE 
2 1 0  IT  =  1 
220  READ  (l;  100)  T 

IF  ( T ( 1  )  +  T ( 3  )  )  230,  250,  230 
230  Dc  2*0  I  =  1,  7,  2 

UEUTD( IT.KSTN)  =  T  < I ) 

F [)(  IT.KSTN)  =  T(  1*1) 

IT  =  IT  ♦  1 
240  CONTINUE 
Go  TO  220 

C  ...  READ  TIRE  SEGMENTS  AND  SPRING  CONSTANTS 
250  IT  =  1 
260  READ  (l;  100)  T 

IF  ( T ( 1 )  ♦  T ( 3  )  )  270,  290,  270 
270  Do  280  I  =  1,  7,  2 

DSEG( IT.KSTN)  =  T(  I) 
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SeGKI I T , KSTN  )  s  T( 1*1) 

it  =  it  ♦  i 
2eu  continue 
Go  TO  260 

c  ...  read  tire  damping  force  profile 

290  It  =  1 

>500  READ  (1;  ICO)  T 

If  ( T ( 1 )  ♦  T ( 3  )  )  310,  150,  310 
310  Do  320  I  =  1,  7,  2 
DEL  D ( IT.KSTN)  =  T ( I > 

FfD<  IT.kSTM  =  T(  1*1) 

it  =  it  ♦  i 
320  continue 

Go  TO  300 

330  continue 

c  ...convert  from  degrees  to  radians 
Do  360  jSj  =  1,  NS US 

c  ...  test  for  sameness 

If  (KSAME(JSJ))  360,  340,  360 
3 4 U  JS  =  JSJ 

NnN  =  NSEG(JSJ) 

Do  35C  I  =  1,  NNN 

DsEGI I , JS)  =  DSEGf I , JS)  •  DTr 

Cxd,JS)  =  S  i  N  !  DSEG  (  1  ,  JS  )  )  •  TREFDUSJ) 

C  y  (  I  ,  JS  )  =  -  COS  <  DSEG  (  I  ,  JS  >  )  •  TREFDUSJ) 

350  continue 
360  continue 

c  ...  read  terrain  profile  table 

NTERK  s  1 

It  =  l 

370  RFAD  (l;  100)  T 

If  (T(1)  ♦  T ( 3 )  )  380.  400,  380 
380  DO  390  I  =  1,  7,  2 
STA< IT)  =  T( I ) 

HgT< IT)  =  T  <  1  *  1  ) 

IT  =  IT  ♦  1 

390  continue 

Go  TO  370 

400  GO  TO  (410,  420),  NTERR 
410  NTERK  S  2 

iterl  =  IT  -  1 
Go  TO  370 

420  ItERR  :  It  -  1 
AlPHAO  =  ALPHA 
A|_PhAX  =  ALPHA  •  DTR 
S 1 N  A  =  SIN(ALPHAX) 

CoSA  =  COS(ALPHAX) 

C  ...  ROTATE  TERRAIN  THROUGH  SLOPE  ANGLE  ALPHA 

DO  440  i  =  1,  ITERR 
IF  (STA ( I ) )  440,  440,  430 
430  continue 

XX  =  STAU)  •  COSA  -  HGT(I)  •  SINA 
YY  =  STA(I)  •  SINA  ♦  HGT(I)  •  COSA 
S  T  A ( I  )  I  XX 
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HOT ( I  )  :  YY 

<40  continue 
betaop  :  beta 
uetax  s  beta  *  dtr 

S ! NB  =  SlN(BETAX) 

c  convert  terrain  profile  to  travel  path  angle  beta 

DO  4 1> 0  1  =  1,  I  T err 
S  T  a ( I  )  :  s  T  A (  I)  /  SlNB 

<50  continue 

c  compute  lag  distances  of  wheel  CENTERS 

T  T  T  o  =  -  L  T  I  R  C 1  )  ♦  L  T  I  R  <  2  ) 

IF  (LOS(BETAX))  <90,  <60,  <Bo 
<60  Do  <70  |  =  1,  NSuS 
<70  WlAGI I )  s  0.0 
GO  TO  5?0 
< 8 U  HlAG(I)  =  O.o 

Wl AG ( 2  )  =  TTTD  /  (SIN(BETAX)  /  COS  <  BE  T AX  ) ) 

Go  TO  500 
<90  W|_AG  (  2  )  =  0.0 

w  l  a  r,  <  i  >  s  tttd  /  <  -  s  i  n  (  be  tax  >  /  cos<betax>) 

50U  K  :  NSUS  -  1 

Uo  510  1  =  3,  K .  2 
J  8  1  I  ♦  1  >  /  2 

WlAGU)  =  (LAX(l)  -  LAX(J))  »  WLAGU) 

WL*G(I*1)  =  f  L  A  X  C 1  )  -  L  A  X ( J  )  )  ♦  WL  AG ( 2 ) 

•jio  continue 
520  continue 


equilibrium  computation 


TEST  FOR  EQUILIBRIUM  COMPUTATION 
IF  (bCGMT)  550,  530,  550 
c  ...  read  initial  values  at  equilibrium 

530  READ  (li  1UO )  C V ( I  )  ,  I  =  1,  N1NVAR) 
do  5<n  i  s  i,  ninvar 
vl  m  s  v  ( i ) 
u<  I  )  =  0 , 

5 4 u  continue 

Go  TO  590 

c  ...  SET  UP  INITIAL  values  FOR  EQUILIBRIUM  COMPUTATION 

55 U  Do  560  I  =  l,  NINVAR 
V  (  1  )  =  VL< I  )  =  0. 

560  Dfl)  =  0. 

D(NINVAR)  =  1.0 

XB  =  bcght 
DO  5/0  1  =  1,  J A X 
570  X A (  I  )  =  ACGHT 
C  ANALYSIS 

C  ...  TOMPuTE  INITIAL  VALUES  OF  ACCELERATION 

call  derIv 

C 

c  ...the  initial  time  step  for  equilibrium  computaions 
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IS  AN  INPUT  VALUE,  THE  MAXIMUM  ALLOWED  IS 
8  TIMES  THE  INITIAL.  THE  MINIMUM  IS  1/8  THE  INITIAL. 

Hh  *  DEL T I M 
HmAX  =  6,  •  HH 
Hh I N  =  .125  •  HH 

...  INTEGRATE 

580  continue 
call  ntgRat 

...  TEST  FOR  FINISH  (EQUILIBRIUM  IS  ASSUMED  TO  BE 
PEACHED  AFTER  SOME  ARBIRARY  time. the  TIME  CHOSEN  M*Y  BE 
CHANGED  AS  APPROPRIATE  for  the  VEHICLE. 

IF  (Time  -  3.  )  5S0,  590.  590 
...  PRINT  EQUILIBRIUM  CONDITIONS 

590  continue 

print  120,  (v( n ,  1  =  1.  ninvari 
. . .  test  f or  continuation 
Print,  "Stop  TImE.STOP  DISTANCE",  ♦  • 

Read*  tstop.  hstop 

IF  (TSTOP. EO.O. )  CALL  EXIT 
MOTION  s  1 


DYNAMIC  RESPONSE  COMPUTATION 


...  SET  UP  INITIAL  VALUES 
V(NINVAR)  s  0.0 
D(NINVAR)  *  1.0 
HP  1  0.0 
HD  *  VEL 
Vl<47)  s  HD 

...  COMPUTE  THE  INITIAL. MAXIHUM,  AND  MINIMUM  TIME 
STEPS  FOR  DYNAMIC  RESPONSE  AS  A  FUNCTION  OF  THE 
FORWARD  VELOCITY  OF  THE  VEHICLE.  THE  MAXIMUM 
ALLOWABLE  FORWARD  DISPLACEMENT  PER  TIME  STEP  WAS 
CHOSEN  TO  BE  8  INCHES, THE  MINIMUM  1/2  INCH, 

DO  600  I  s  1,  10 

Hh  «  3 ,  •  2 .  ••(-!) 

IF  (HH  -  8.  /  VEL)  610,  600.  600 

6oo  continue 

610  HmAX  =  HH 

HmIN  S  HMAX  /  16. 

WRITE  (3)  ICOMNT,  VEL.  NSEG(l).  HMAX,  H«IN 

call  output 

...now  ready  to  simulate  forward  motion  of 

the  VEHICLE  OVER  A  TERRAIN  OR  OBSTACLE  COURSE. 

620  call  ntgrat 
call  output 

C  ...  TEST  FOR  FINISH 
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ir  (Time  -  tstop>  63o.  640,  640 
630  If  (HP  -  MSTOP)  620,  6-10,  640 
64 U  CONTINUE 

print,  "Run  completed" 
call  exit 
end 
c 

c  . . .  subroutine  ntgrat 

...this  subroutine  uses  runge-kuuta-mersons 
integration  forulae.  the  step  size  is  automatically 
adjusted  between  the  cohputed  limits  based  on  the 
computed  integration  ERROR  FOR  The  VERTICAL 
motion  Of  THE  AXLES.  THE  STEP  SIZE  IS  ALSO  REDUCED 
IF  THE  DFRIVITIVE  ROUTINE  RETURNS  An  INDICATION 
THAT  DISPLACEMENTS  or  velocities  exceed  THE  MAXIMUM 
values  provided  as  vehicle  PARAMETERS, 

Subroutine  ntgrat 
iou  luERR  =  i 

HH03  :  hh  /  j. 

Do  UO  I  =  1 ,  N  1  N  V  A  R 
_  Pi ( I )  s  Hh03  a  D ( I ) 

1 1  u'  V  (  I  >  =  VL  (  I  )  ♦  Pi  (  I  ) 

call  derIv 

c  ...  luERR  is  returned  as  2  if  a  table  look-up  error  OCCURS 
Go  to  <121,  260),  LUERR 

120  continue 

00  130  I  :  1,  NINVAR 
P2< I  >  *  HH03  «  D ( I  ) 

1 3 U  V ( 1 )  =  VL  (  1 )  ♦  (PI (  I  )  ♦  P2( I  )  )  •  .5 

call  derIv 

go  to  (140,  280),  LUERR 

140  continue 

do  150  i  =  l.  ninvar 
P 3 ( I  )  =  HH03  o  D (  I  ) 

150  V(l)  s  VL  (  I  )  ♦  Pl(I)  •  .375  ♦  P3 (  I  )  *  1.125 

call  deriv 

GO  TU  (160,  280),  LUERR 

i6 u  continue 

do  i 7 u  i  =  i,  ninvar 
P 4  < I  )  =  HH03  «  D (  I  ) 

170  V(l>  =  VL ( I )  ♦  P 1 (  I  )  •  1.5  -  P  3  <  I  )  •  4,5  ♦  P4  (  I  )  o  6.0 

call  deriv 

Go  to  ( 1 6 0  ,  2  8  0)  ,  LUERR 
ieo  continue 

Do  1V0  I  s  1,  NINVAR 
P 5  <  I  7  =  Hh03  a  D (  I  ) 

Du  =  (PI  ( I  )  ♦  P4 ( I  )  •  4 .  ♦  P5(  I  >  )  •  .5 
19 U  V ( I  )  :  VL  (  1  )  ♦  Du 
lDOUbL  =  1 
Do  220  1  =  1,  J  A  x 

K  =  1  ♦  6 

ERR  =  ADS((Pl(K)  -  P  3  <  K  )  •  4,5  ♦  P  4 ( K  >  • 
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ooooooo  no  o  noooo  ociooooo* 


4.0  -  P 5 ( K )  *  .5)  *  .2) 

...bRROR  OF  10  IN/SEC  IS  THE  LIMIT  BEYOND  WHICH  A  STEP 
S I  Zb  REDUCTION  OCCURS. 

...bRRoR  OF  1  In/SEC  IS  THE  LIMIT  BELOW  WHICH  A  STEP 

s i zb  increase  occurs. 

if  (Err  -  10 . >  200.  200.  2flo 
2ou  ir  ( bRR  -  1.)  2 20.  210.  210 
210  iDOUbL  -  0 
220  continue 

IF  (IOOUBU  250.  250  ,  230 
230  Ir  (HH  -  HMAX)  240,  250.  250 
240  Hu  =  HH  *  2. 

250  NxTIME  =  5 
CALL  DERlV 

Go  TU  <26u.  280).  LUERR 
260  CONTINUE 

Do  2?0  I  =  1,  NINVAR 
Vl<  I  )  s  V ( 1 ) 

2 7 u  continue 
return 

280  Ir  (hH  -  hMIN)  300.  300.  290 
290  Hh  =  HH  •  .5 

Go  TU  100 

300  Print,  "STEP  SIZE  LESS  THAN  MINIMUM, RUN  TERMINATED" 
call  EXIT 
End 

. . .  SUBROUTINE  DERIV 

...  this  SUBROUTINE  IS  USED  BY  NTGRaT  TO  COMPUTE  DERIVlTlVES 
subroutine  deriv 

IF  (MOTION)  100.  100,  450 

...  MOTION  is  0  FOR  EQUILIBRIUM  COMPUTATION  ,  1  FOR 
dynamic  computation. 


DERIVATIVE  ROUTINE  FOR  EQUILIBRIUM  COMPUTATION 


IF  (MOTION)  100,  100,  450 

loo  continue 

C  ...  CALCULATE  DEFLECTIONS  OF  TIRES  AND  SPRINGS 

Uo  IPO  jSj  =  1,  NSUS 
J  A  J  *  (JSj  ♦  1)  /  2 
C  ...  TIRE  DEFLECTIONS 

DL  T ( JS  J )  =  TREFD(JSJ)  -  (XA(JAJ)  -  LTIR(jSj)  •  S I N ( PH A ( J A J ) ) ) 
IF  (ULT(JSJ)J  110.  120,  120 
110  DlT(JRJ)  =  0.0 


\ 

t 


Cl  po  C» 


CONTINUE 

,  .  ,  SPRING  DEFLECTIONS 
DlS(USJ)  =  SREFD(JSJ)  -  -  UX(JAJ  • 

SIN(THH)  -  LSUS(JSJ)  <*  S I N { G0  >  -  X  A  <  J  A  J  ) 
LASUS(JSJ)  *  S I N ( PHA (  J  A  J  )  )  ) 

...  calculate  deflection  velocities 
IF  (ULT(JSJ))  130.  130.  HO 
DlIDUSJ)  =  o.o 
Go  TO  150 

CONTINUE  .  .  rnc-D 

DlTDIJSJ)  =  -  XDA(JAJ)  *  LTIR(JSJ>  •  COovP 

HA ( JAJ ) )  0  PH0A ( JAJ ) 

CONTINUE  _  rnCITHBI  • 


DLSDUSJ)  -  '  XDB  ♦  L  A  X  (  J  A  J  )  *  COSCTHBJ  •  _ 

THUB  ♦  LSUS(JSJ)  •  COS(GB)  •  GDB  *  *D*tJ*J) 
LASUS(JSJ)  0  COS(PHA( JAJ) )  «  PHD  A ( J A  J ) 

CONTINUE 

DO  400  JSJ  =  1»  NSUS 
JAJ  =  (JSJ  ♦  1)  /  2 
...  TEST  f CR  SAMENESS 
IF  (KSAMfc(JSJ))  170.  180.  170 
JS  =  KSAME(JSJ) 

Go  TO  190 
JS  =  JSJ 
continue 

N^,  look  UP  SUSPENSION  SPRING  FORCE 

??  DLS.JSJ))  ?00.  20°.  220 

CONTINUE 

fSJmt2;";  Of  "  TLU  ERROR  -  SUSPENSION  SPRING-U, 
LuERH  =  2 

SPSUSJ)  =  FS(l-l.JS)  ♦  (DLSIJSJ)  -  DILJAU-I  _ 

,  JS )  )  o  (FS(I.JS)  -  FS(l-l»JS)>  /  IDEL'A(l.JS) 

DEL  T  A  (  I-l.JS)  )  rnDrc 

...  look  UP  SUSPENSION  DAMPING  FORCE 

??  ruELioilOSI4-  OLSOUSJ,,  23°.  230.  250 
CONT 1 NUE 

E0RMAT2(0/  5X^  "  TLU  ERROR  -  SUSPENSION  DAMP" 1 3 ) 
LUfcRR  =  2 

-  EDM-l.JS,  *  ."LSD.JSJ.  -  . 

,  JS ) >  .  (FD(I.JS)  -  F  D ( I  — 1 »  JS ) >  /  COELTD(I.JS) 

DEL  T  D ( I-l.JS) )  „  „ 

...  COMPUTE  TIRE  SPRING  FORCE 
SpT(JSJ)  =  0. 

IF  (ULT(JSJ))  360.  360.  260 
j  continue 

CC  c  2 !  =  XA(  JA J > S- ’ L  T I R ( JSJ)  •  S I N ( PH A < JA J > > 
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LOCATE  ALL  POINTS  ALONG  TIRE  CIRCUMFERENCE 
Do  270  I  =  1,  NNN 
RXY(1,  I  )  =  XI  ♦  CXI  I , JS) 

RX Y  I  2  ,  I  )  =  Y 1  ♦  CYI! , JS) 

27U  CONTINUE 

, . .COMPUTE  DEFLECTIONS 
Do  34(1  I  =  1,  NNN 
Sh  =  CYI I , JS)  /  CXI  I . JS) 

IT  (Sm)  290,  280,  290 
28U  XX  =  1  .  E3t> 

Go  TO  300 

290  XX  =  I  -  Y1  ♦  SM  •  XI)  /  SH 
JOU  Yy  s  0, 

X2  =  0, 

X3  =  2.  •  TREFD(JSJ) 

KkERK  =  KBTWN(X2,XX,X3) 

GO  TO  (310,  340),  KKERR 
3 1 U  KkERT  =  KB  TWN  <  XI , XX, RXY 11,1)) 

GO  TO  (320,  340),  KKERT 
320  KkERU  =  KBTWN(Y1,YY,RXYI2,I)> 

Go  TO  (330,  340),  KKERU 
330  Rx Y  ( 1 , I  )  =  XX 
RXY(2,  I  )  =  YY 

340  Continue 

Do  3t> 0  I  =  1,  NNN 
A  A ( 1 >  =  CC(l)  -  RXY  <1,1  ) 

A A  I  2 )  =  CCI2)  -  RXY I  2 ,  I  ) 

T  M *G I  I  )  =  SQHTIAA(l)  «  AA(1)  ♦  A A  C  2  >  •  AAI2)) 

RX Y ( 2 ,  I  )  =  AA ( 2 )  /  TMAGI 1  ) 

TmAGII)  =  (TREFO(JSJ)  -  TMAG(I))  •  SEGKH.JS) 
SpT(JSJ)  =  SPTIJSJ)  ♦  TMAGII)  •  RXY  1 2 . I  ) 

350  CONTINUE 

360  continue 

...  LOOK  UP  TIRE  DAMPING  FORCE 
Do  370  I  :  1 ,  24 

IF  I DELD I  I , JS )  -  DLTDIJSJ))  370,  370,  390 
370  CONTINUE 

PRINT  380,  JSJ 

380  format  <  /  5x>  "tlu  error  -  TIRE  damp»i3) 

LuERR  =  2 
RETURN 

390  DMT(JSJ)  =  FFD I  I -1 , JS  )  ♦  IDtTDIJSJ)  -  DELD ( I -1 

.  JS  )  )  «  (l-FDII.JS)  -  FFDI I -1 ,  JS  >  >  /  I  DELDI I ,  JS  ) 
DELDI I - 1 »  JS ) ) 

. . .  COMPUTE  TOTAL  FORCES 
F(JSJ)  =  SPS(JSJ)  ♦  OMSIJSJ) 

FFUSj)  =  SPTIJSJ)  ♦  DMTIJSJ) 

FlATIJSJ)  =  F(JSJ)  *  SIN(GB) 

4oo  continue 

COMPUTE  SUM  of  SUSPENSION  FORCES 
SuMB  =  0.0 

Do  410  JSJ  :  1,  NSUS 


■H.J"  ‘  nuw»»n.  )  'U.ii,i  ,J|  W 
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*1U  SjjMB  s  SUM6  ♦  F(JSJ) 

c  ...  compute  body  acceleration  and  velocity 

D  { 1 )  =  (SUHB  -  WTB)  <*  GG  /  WJB 
0(2)  =  v  < 1 ) 

c  ...  compute  sum  of  body  pitch  MOMENTS 
SuMOM  =  o.c 
Do  420  JSJ  =  1,  NSUS 
JAJ  (jSj  ♦  1)  /  2 

420  SuMOM  =  SUMOM  -  LAX(JAJ)  *  F{JSJ)  •  COS(THB) 

C  ...  COMPUTE  HODy  pitch  angular  ACCELERATION  AND  VELOCITY 
D(J)  =  SUMOM  /  IBP 
D  (  4  )  =  v  (  3  ) 

C  ...  COMPUTE  SUM  OF  BODY  ROLL  MOMENTS 

SUMOM  =  0.0 
Do  43o  JSJ  =  1,  NSUS 

430  SuMOM  =  SUMOM  -  LSUS(JSJ)  «  CO  S  <  G  B )  •  F(JSJ) 

C  ...  COMPUTE  BODY  ROLL  ANGULAR  ACCELERATION  AND  VELOCITY 

Dp)  =  SUMOM  /  IBR 

D ( 6  )  =  v< 5) 

C  ...  COMPUTE  AXLE  ACCELERATION  AND  VELOCITY 

Do  4 A o  jSj  =  1,  NSUS,  2 
JAJ  s  (jSj  ♦  1)  /  2 
C  DOH 

D  (  J  A  J  ♦  6  )  =  (FF(JSJ)  ♦  FFUSJ  +  l)  -  F<jSj)  - 
&  F ( JSJ*1 )  -  W  T  A  x ( J  A  J  )  )  #  GG  /  WTAX(JAj) 

D ( J A J ♦ 16  )  =  V  (  J  A  J ♦ 6 ) 

C  ...  COMPUTE  AXLE  ROLL  ANGLE  ACCELERATION  AND  VELOCITY 

D ( JA J*26 )  =  (  -  FF(JSJ)  *  LTIR(JSJ)  *  COS(P 
&  HA(JAJ)>  -  FF(JSJ*i)  •  LT I R ( JS  J*1 >  °  COStPH 

4  A(JAJ))  *  F(JSJ)  <•  L  ASUS  (  JS  J )  •  COStPHAI  JAJ)  )  ♦ 

4  F ( JS J ♦ 1 )  •  L ASUS  <  JSJ*1 )  •  COStPHAt JAJ)  ) >  /  IAXR(JaJ) 

D ( J A J  +  36  )  s  V ( J 4  J  +  26  ) 

440  Continue 

D  (  4  7  )  =  0.0 
D ( 4  8 )  =  0.0 
RETuHn 

C«»»o»o*.  DERIVATIVE  ROUTINE  FOR  DYNAMIC  RESPONSE  COMPUTATION  • 

C  O«O«oOo«OOO»ftoe«O«»O«6»«O|«*0O»«t««IM«4)«eiO»««O|MOH*Ot 

450  CONTINUE 
LuERR  =  1 
FhOrIz  =  0.0 
HoRmOm  =  0.0 

C  ...  COMPUTE  DEFLECTIONS  OF  TIRES  and  TIRE  FORCES 

Do  620  JJJ  =  1,  NSUS 
JSJ  =  NSUS  -  JJJ  *  1 

c  . . .test  for  sameness 

If  (KSAMEiJSJ))  460.  470.  46o 
460  JS  =  KSAME(JSJ) 

Go  TU  480 
4 7 U  JS  =  JSJ 

480  continue 

c 

JAJ  =  (JSJ  +  u  /  2 


FFV(  1 )  =  0.0 
F  F V ( 2  )  =  0.0 
NNN  *  NStG(JSJ) 

V V V  < 1 )  =  0. 

VvV( 2)  =  0. 

WRAD  =  TREFD(JSJ) 

...  FIND  EXTREME  REAR  OF  TIRE  AND  EXTREME  FRONT  OF  TIRE 
XwCN  =  HP  -  WLAG(JSJ) 

XRET  =  XWCN  -  HR AD 
XFET  =  XKCN  +  WRAD 
...  LOOK  UP  TERRAIN  COORDINATES 
KLR  =  JSJ  /  J A J 
GO  TO  ( 4 9 o .  500  ),  KLR 
KtRI  =  1 
K  tR2  =  ITERL 
Go  TO  510 
KtRI  =  ITERL  ♦  1 
KTR2  =  ITEHR 

continue 

Do  520  I T R  =  KTRl,  K TR2 

IF  (XRET  -  STA(ITR))  530,  520,  520 

continue 

print,  "vehicle  is  seyond  limits  of  the" 

,  "  TERRAIN, RUN  TERMINATED" 

call  exit 
CONTINUE 
call  segwhl 

IF  (LUERR.E0.2)  RETURN 

...  COMPUTE  WHEEL  center  velocity  vector 

VEL V ( 1 )  =  HD 

VELV12)  =  XOA(JAJ)  -  LTIR(JSJ)  •  COS(PHA(JAJ))  •  PHdA(JAJ) 
...  COMPUTE  TIRE  DEFLECTION  VELOCITY 
SDLD  =  -  (  VEL  V  (  1 )  •  VVV  <  1 )  +  VEL  V  (  2  )  •  WV(2)) 

.  .  .  TEST  F OR  SAMENESS 
IF  (KSAMb(JSJ))  540.  550,  540 
JS  =  KSAME(JSJ) 

Go  TO  560 
US  -  JSJ 
CoNT InUE 

...  LOOK  UP  TIRE  DAMPING  FORCE 
Do  5/n  I  =  1,  24 

IF  (DELD(I'JS)  -  SDLD)  570,  570,  600 

conti nue 

continue 

print  590,  jsj,  sold 

luerh  =  2 

return 

Format  C  /  5x<  "TLU-ERROR  -  TIRE  CAMP"  J  3 ,  2X ,  "SDLD*''F10 .3) 
IF  (l  -  1  )  580  ,  500  ,  610 

TDFOK  =  FFD(I-l.JS)  ♦  (SDLD  -  DELD ( I -1 , JS  ) )  * 

(  FF  0  (  1  ,  J  S  )  -  FFD(I-IiJS))  /  (OELDU-JS)  -  DELD  (  I -1 ,  JS  )  ) 
...  COMPUTE  FORCE  VECTOR  AND  ACCUMULATE 
SuPFUR  =  TSFOR  *  TDFOR 
FoRVfcC(l)  =  SUMFOR  •  VVV(l) 


. . . 
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F  0«VEC ( 2 )  =  SUMFOR  *  VV V { 2 ) 

F  F V  < 1 )  =  FORVEC(I)  +  FFV(1> 

F  F V  <  2  )  =  FORVEC ( 2  )  ♦  FFV(2> 

C  ...  SET  VERTICAL  force  component 
FF(JSJ)  =  F  F  v ( 2  ) 

c  ...  accumulate  horizontal  force  component 

F  HOR 1 Z  =  FHORIZ  ♦  FFv(l) 
c  ...  accumulate  horizontal  moments 

HoRMOM  =  HORMOM  ♦  FF V ( 1 )  •  (XB  -  (XA(JAJ)  - 
4  SIN(PHA( JAJ) )  •  LTIR(JSJ))) 

62 u  continue 

c  ...  COMPUTE  DEFLECTIONS  of  springs  and  spring  forces 
Do  760  jSj  -  1.  NSUS 
JAJ  e  (JSJ  ♦  1)  /  2 
IF  (KSAMEtJSJ))  630.  640.  630 
630  JS  =  KSAME(JSJ) 

Go  TO  650 
640  JS  =  JSJ 

650  continue 

c  ...  calculate  spring  deflections 

DlS(JSJ)  s  SSEFD(JSJ)  -  <  X8  -  LAX(JAJ)  • 

4  SIN(ThB)  -  lSUS(JSJ)  •  SIN(GB)  -  XAtJAj)  ♦ 

£  LASUS(JSJ)  •  SIN(PHA( JAJ) ) ) 

C  ...  CALCULATE  the  DEFLECTION  nelocities 

DlSD(JSJ)  =  -  XDH  ♦  LAX ( JAJ )  «  COS ( T  HB  )  • 

4  THDB  ♦  LSUS(JSJ)  o  COS(GB)  •  GD8  *  XDA(JAJ)  - 

4  LASUS(JSJ)  »  COS(PHA( JAJ) )  ®  PHDA(JAJ) 

C  ...  LOOK  UP  SUSPENSION  SPRING  FORCE 

Do  600  I  s  1,  24 

IF  ( DEL T A ( I . JS )  -  DLS ( JSJ ) )  660.  660,  690 

660  continue 
670  continue 

print  680,  jsj 
LuERR  =  2 
return 

680  FORMAT  <  /  5X.  "TLU  ERROR  -  SUSPENSION  SPRING"I3) 

690  if  (I  -  1)  670,  670.  700 

700  SpS(JSJ)  :  FS(I-I.JS)  ♦  <  DL  S ( JSJ  )  -  DEL  T  A  t I -1 

, JS )  )  »  (FS(I.JS)  *  FS(I-l.JS))  7  (DELTA!!, JS)  - 
DELTA!  1-1, JS)  ) 

...  LOOK  UP  SUSPENSION  DAMPING  FORCE 
Do  710  I  =  1,  24 

IF  (DELTD(I.JS)  -  DLSD(JSJ))  710,  710,  740 

7io  continue 
720  continuf 

print  730,  jsj 
LUERR  s  2 
RETURN 

/30  format  <  /  5x.  "tlu  error  -  suspension  damp"i3) 

/ 4  0  IF  (I  -  1  )  720  .  720  ,  750 

750  DMS(JSJ)  =  FD(I-l.JS)  ♦  (DLSD(JSJ)  -  DEL  T  D  < I -1 
4  ,JS>>  «  (fD(I.JS)  '  FD(I-1,JS)>  /  (DELTDCI.JS)  - 

4  DEL TD < I -1 . JS )  ) 

C  ...  COMPUTE  TOTAL  FORCES 


! 

I 

I 

F(JSJ)  =  SPS(JSJ)  +  DMS(JSJ) 

FvERtJSJ)  =  F(JSJ)  •  COS(GB) 

FLAT(JSJ)  =  F(JSJ)  •  SIN(GB) 

/6U  continue 

c  compute  sum  of  suspension  forces 

Sums  =  o.u 

Do  770  JSJ  =  1,  NSUS 
770  SUMB  =  SUMB  ♦  FVER(JSJ) 

C  COMPl  te  body  acceleration  and  velocity 

D ( 1 )  =  (SUMB  -  WTB )  *  GG  /  WTB 
D(2)  =  XDB 

c  ...  compute  sum  of  body  PITCH  moments 

SuMOM  a  0.0 
Do  700  JSJ  =  1.  NSUS 
JAJ  =  (JSJ  ♦  1)  /  2 

780  SuMOM  =  SUMOM  -  LAX(JAJ)  •  F(JSJ)  #  COS(THB) 

C  PALO 

SuMOM  a  SUMOM  -  HORMOM 

C  ...  COMPUTE  BODY  PITCH  ANGULAR  ACCELERATION  AND  VELOCITY 
D(J)  s  SUMOM  /  IBP 
D ( A )  =  THDB 

C  ...  COMPUTE  SUM  OF  BODY  ROLL  MOMENTS 

SuHOM  =  0.0 
DO  7 V 0  JSJ  =  1.  NSUS 

790  SuMOM  =  SUMOM  -  LSUS(JSJ)  •  F(JSJ)  •  COS(GB) 

C  ...  COMPUTE  BODY  ROLL  ANGULAR  ACCELERATION  AND  VELOCITY 

D  ( t»  >  =  SUMOM  /  IBR 
0(6)  =  GOB 

c  ...  COMPUTE  AXLE  ACCELERATION  AND  VELOCITY 

Do  800  JSJ  =  1,  NSUS.  2 
JAJ  =  (JSJ  ♦  1)  /  2 

D ( JA  J ♦ 6  )  r  (FF(JSJ)  ♦  FF ( J  S  J ♦ 1 )  -  F(JSJ)  - 
4  F ( JS J ♦  1 )  -  WTAX(JAJ))  *  GG  /  WTAX(JAJ) 

D ( JAJ  +  16  )  =  XDA(JAJ) 

C  ...  COMPUTE  axle  roll  ANGLE  ACCELERATION  AND  VELOCITY 
D( JAJ+26)  =  (  -  FF (JSJ I  •  LTJR(JSJ)  *  COSIP 
«  HA(JAJ))  -  FF(JSJ*1)  *  L  T  I R { JS J ♦ 1 )  •  COS(PH 

4  A(JAJ))  *  F (JSJ)  •  LASUS(JSJ)  •  COS ( PHA ( JA J )  >  ♦ 

4  F ( JS  J ♦ 1 )  •  L ASUS  <  JS J* 1 )  •  COS  <  PHA ( JA J  )  )  )  /  IAXR(JaJ) 

D(JAJ  +  36  )  =  PHDA(JAJ) 

eoo  continue 

C  ...  COMPUTE  HORIZONTAL  ACCELERATION  OF  TRUCK 

D(47)  s  ACCEL 
D ; A 8 )  s  HU 
C  ...  EXIT 

return 
End 
c 
c 

c  . . . .  SUBROUTINE  OUTPUT 

c 
c 

c  ...  THIS  SUBROUTINE  WRITES  ACCELERAT I ONS , VELOC I T I ES . 

C  AND  DISPLACEMENTS  TO  AN  OUTPUT  FILE  AT  EACH  INTEGRATION 
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c  sTEp.IT  is  CALLED  AFtEr  EACH  $UCES5FUL  INTEGRATION  STEP. 

C  IT  CAN  BE  CHANGED  TO  OUTPUT  THESE  VARIABLES  IN  ANY 
C  DESIRABLE  FORM.  ANY  OTHER  VARIABLES  OF  INTEREST  MAY 

C  BE  OUTPUT  HFRE  SINCE  ALL  SIGNIFICANT  VARIABLES  ARE 

c  in  Common. 

c 

subroutine  output 

1 NDX  =  I NDX  ♦  1 
C  ...  IS  THE  FILE  FULL 

IF  <INDX,lT.i26)  GO  TO  120 
c  . . .yes. create  a  new  file 

call  closeftj) 

NmF  L l 2 )  =  NMFL ( 2  )  ♦  1 
lOU  FoHMmT  (A6) 

print  ho.  pltout.  vi5d 
no  format  (lx.  A6,  ix.  "begins  at",  ix.  F9.5> 
write  (3)  DUMMY 
call  closef is.  pltout) 

Call  OPENFI3,  PLTOUT) 

C  ...  WRITE  THE  HEADING  INFO  AT  the  BEGINNING  OF  Each  nt-E, 
WRITE  (3)  ICOMNT,  SPD.  NSEG ( 1 )  i  MMAX,  HMIN 
InDX  =  1 

120  write  (3)  D.  V 
RF  T  URN 
EnU 
c 
c 

c  «ooo»ooo#ooooooo«»oo  function  kbtwn 

c 

c  ...  this  function  is  USED  By  SUBROUTINE  segwhl  TO 

C  DETERMINE  IF  a  POINT  IS  BETwEEN  two  other  points. 

c 

FuNCT ION  kBTWNIA,  P,  B) 

IF  (A  -  P)  110,  120,  100 
100  IF  (b  -  P)  120.  120.  130 
c  IS  P  BFTwEEn  A  and  b 

HO  IF  (B  -  P)  130.  120.  120 
C  VALID  return 

120  KRTWN  =  1 

re  tuhn 

C  INVALID  RETURN 

130  K B  T  w iv  =  2 
ReTuKn 
End 
c 
c 

C  iimmiiniimmii  SUBROUTINE  SEGWHL 

c 

c  ...this  subroutine  is  usfd  by  the  dynamic  response 

C  PORTION  OF  SUBROUTINE  DERIV  TO  COMPUTE  TIRE  DEFLECTION 

C  FORCE  VECTORS  USING  THE  SEGMENTED  TIRE  CONCEPT. 

C 

subroutine  segwhl 

cell)  =  XwCN  -  STA  <  I  TR-1  ) 


r 


C C ( 2 )  =  XA(JAJ)  -  LTIR(JSJ)  *  S I N ( PHA ( J A J > ) 

XI  =  CC(1) 

Y 1  =  CC ( 2 ) 

C  ...LOCATE  all  points  ALONG  TIRE  CIRCUMFERENCE 
Do  1UO  1  =  1,  NNN 
RX Y  C 1 , I )  =  XI  ♦  CXI  I  , JS) 

RxY(2, I )  =  Y 1  ♦  CYI  I  .JS) 

lot)  continue 

c  ...Intersect  segment  centers  witr  terrain 

MOLE  =  1 

DO  220  I  =  1.  NNN 
DmIN  =  1.E10 

Sm  =  CYI 1 , JS)  /  CXI  I  . JS) 

KTR  =  ITR 
X2  =  0. 

Y2  =  HGTIKTR-1) 

liu  continue 

X3  =  STA(KTR)  -  STAIITR-l) 

Y  3  =  HGTIKTR) 

ST  =  I Y2  ••  Y3)  /  (X7  -  X3> 

IF  (SM  -  ST>  130,  120.  130 
12U  Xx  =  1.E35 
GO  TU  140 

130  Xx  =  I Y 2  -  Yl  ♦  SM  «  XI  -  ST  •  X2)  /  I SM  -  ST) 

140  Yy  =  ST  *  (XX  "  X2)  ♦  Y2 
KERR  =  kBTRN(X2.XX.X3> 

GO  TO  ( l50 1  180  ),  KERR 
150  KERS  =  KBTKNIY2.YY.Y3) 

GO  TO  (160,  180),  KERS 
160  KERT  j  KBTWNIXl , XX.RXYI 1, I ) > 

GO  TO  I  1 7  o .  180),  KERT 
170  KERlJ  =  KB  TWN  I  Yl .  YY.RXYI2,  I  )  ) 

GO  TO  (200,  180),  KERU 
180  IF  (XFET  -  STAIKTR))  220,  220,  190 
190  KTR  =  KTR  *  1 
X2  =  X3 

Y  2  =  Y3 
GO  TO  HO 

200  DmINI  =  SORT(IXX-Xl)  •  I  XX  -  Xl>  ♦  <YY  -  Yl)  •  (YY  -  YlM 
IF  (bMlNl  -  OMlN)  210.  210.  180 
2 1 0  DmIN  —  DMIM 
RXY I  1 .  I )  =  XX 
RX Y I  2  .  I  )  r  YY 
MOLE  =  0 
Go  TU  180 

22U  continue 

IF  (MOLE)  290.  290.  230 
230  Continue 

KTR  =  ITR 

240  If  (STA(KTR)  -  XWCN )  250.  260.  260 
250  KtR  *  KTR  ♦  1 
Go  TO  240 

260  tempi  =  hgtiktr-d 

T  E  MP  2  =  STA(KTR-l) 
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PRHUB  :  TEMPI  ♦  ((HGT(KTR)  -  TEMPI)  /  <STA( 

&  kTR)  -  TEMP?))  •  (XwCN  -  TEMP2) 

IF  (VI  -  PBHuB)  270.  270,  290 
27U  print  2b0,  jsj.  xwcn 

280  TofiMAT  (  //  IX.  "CENTER  OF  WHEEL  NUMBER", 

«  13.  IX.  "IS  below  THE  TERRAIN  SURFACE  AT  XWCN  s»,  F 1 2 . 2  > 

LUERh  =  2 
Hr  TURN 

290  continue 

c 

UO  300  I  si,  NNN 

RXY(1. I  )  =  CC( 1 )  -  RXY(1,  I  ) 

RxY(2. I )  =  CC(?)  -  RXY ( 2  .  I  ) 

3oo  Continue 

C  ...NORMALIZE  RXV  ANT)  let  TMAG  BE  MAGNITUDE 
Do  31C  I  si,  NNN 

A  A  (  1  )  s  K  X  V  (  1  ,  I  ) 

aa ( 2)  =  RxY ( 2  .  !  ) 

T  M AG ( I  )  =  SORT  ( A  A  (  1  )  •  A  A (  1  >  ♦  A  A ( 2 )  •  A  A ( 2  >  > 

RxY (1,  I  )  :  AA (  1  )  /  T  M  A  G ( 1  ) 

R  X  Y { 2 ,  I  )  s  AA (2)  /  T  MAG ( I  > 

3 i o  continue 

C  ...COMPUTE  FORCE  FOR  EACH  SEGMENT 

Do  320  I  si.  NNN 

Tmag(I)  =  (hraD  -  TmaG(I))  •  SEGK(l.jS) 

320  CONTINUE 

C  ...LOMFUTE  COMPOSITE  FORCE  VECTOR 

UO  330  I  S  1,  NNN 

WVfl)  :  V  V  V  (  1  )  ♦  TMAG(I)  •  RXY(l.I) 

VvV<<)  r  VVV(?)  ♦  T  H  A  G ( I  >  •  RXY(?,!> 

3 3 u  Continue 

C  ...COMPUTE  FORCE  MAGNITUDE.  TSFOR 

Tsfok  s  SgRT(vvv(i)  »  vwu)  .  vvv  ( 2 )  •  v  v  v  <  2  )  ) 

IF  (TSFoR)  350.  350.  340 

340  continue 

VvV(l)  :  V  V V ( 1  )  /  TSFOR 
VvVU)  =  VV V ( 2 )  /  TSFOR 
GO  TO  360 
350  V  v V  < 1  )  =  u. 

V  V  V  (  2  )  s  1, 

36U  continue 
return 
end 
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